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This  document  (Volume  IV)  is  the  Theoretical  Manual  of  the  Aeroservoelastic  (ASE) 
interaction  module  developed  to  facilitate  ASE  analysis  and  the  application  of  ASE  stability 
and  response  constraints  within  ASTROS. 

At  AFRL/Wright-Patterson,  Captain  Gerald  Andersen  was  the  contract  monitor  and  Dr.  V.  B. 
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and  others  from  AFRL  during  the  course  of  the  present  phase  on  the  development  of 
ASTROS*  are  gratefully  acknowledged. 
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Chapter  1 
Introduction 


The  common  approach  for  formulating  the  dynamic  equations  of  motion  of  aeroelastic  sys¬ 
tems  uses  normal  modes  of  the  structure  as  generalized  coordinates.  Control-surface  de¬ 
flection  modes  may  be  added  when  the  interaction  with  the  control  system  is  considered. 
Complex  gust  velocity  modes  may  also  be  added  to  analyze  the  response  of  the  structural 
and  control  systems  to  continuous  gust.  The  unsteady  aerodynamic  force  coefficients  are 
defined  with  respect  to  these  modes. 

The  unsteady  aerodynamic  codes  in  ASTROS  assume  that  the  structure  oscillates  har¬ 
monically.  Transcendental  unsteady  aerodynamic  matrices  are  calculated  for  various  reduced 
frequency  values.  The  flutter  and  gust-response  modules  in  ASTROS  use  second-order  for¬ 
mulations  for  stability  solutions,  frequency  response,  and  frequency-domain  control  synthesis. 

The  application  of  various  modern  control  design  techniques,  simulations  and  optimiza¬ 
tion  procedures  require  the  aeroservoelastic  equations  of  motion  to  be  transformed  into  a 
first-order,  time-domain  (state-space)  form.  This  transformation  requires  the  aerodynamic 
matrices  to  be  approximated  by  rational  functions  (ratio  of  polynomials)  in  the  Laplace  do¬ 
main.  The  order  of  the  resulting  state-space  model  is  a  function  of  the  number  of  selected 
modes,  the  number  of  aerodynamic  approximation  roots,  and  the  approximation  formula. 
The  main  considerations  in  constructing  the  model  are  its  size  (which  affects  the  efficiency 
of  subsequent  analyses),  its  accuracy,  and  the  model  construction  efforts. 
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Most  commonly  used  rational  aerodynamic  approximation  methods  are  variations  of  ei¬ 
ther  Roger’s  classic  approximation  [1],  which  is  based  on  term-by-term  least-square  fits  with 
common  denominator  roots,  or  the  minimum-state  (MS)  method  of  Keirpel  [2,  3]  which  is 
based  on  a  more  general  approximation  function  with  coupling  between  terms.  Consequently, 
the  MS  procedure  requires  computationally  heavier,  iterative,  nonlinear  least-squaxe  solu¬ 
tions.  Recent  developments  of  the  MS  method  [4-8]  improved  the  accuracy,  flexibility  and 
computational  eflSciency  of  the  MS  scheme.  The  number  of  aerodynamic  states  added  by 
the  MS  method  in  realistic  aeroelastic  design  models  is  typically  6-8  times  smaller  than 
those  added  by  Roger’s  approximation  with  the  same  level  of  model  accuracy  in  subsequent 
analyses.  This  makes  the  MS  method  very  attractive  in  repetitive  aeroelastic  optimization 
studies  which  are  based  on  modes  of  a  baseline  structure,  and  in  cases  which  involve  control 
synthesis. 
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Chapter  2 

GenereJized  Matrices 


The  modal  approach  to  structural  optimization  is  based  on  using  a  set  of  low-frequency 
normal  modes  of  the  baseline  structure  a  fixed  set  of  generalized  coordinates  throughout 
the  entire  optimization  process,  or  at  least  for  a  major  optimization  cycle  composed  of 
several  design  steps.  The  implementation  of  the  modal  approach  in  ASTROS  required  a 
slight  modification  of  this  approach.  Instead  of  keeping  the  modal  coordinates  fixed,  they 
axe  chajiged  in  each  design  step,  but  the  new  set  of  modal  coordinates  is  assumed  to  be  a 
linear  combination  of  the  baseline  set.  Mathematically,  there  is  no  difference  between  the 
two  approaches. 

It  is  currently  assumed  that  the  model  does  not  include  extra  points,  such  that  the  h-set 
and  i-set  are  identical.  If  extra  points  would  exist,  the  formulation  below  would  apply  to  the 
n  blocks  of  the  dynamic  matrices  only,  with  the  other  blocks  remain  unchanged  during  the 
design  process.  We  also  assume  here  that  the  contributions  of  the  design  variables  to  the 
stiffness  and  mass  matrices  are  linear. 
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2.1  The  Modal  Assumption 


The  set  of  calculated  baseline  modes  [^oi]>  defined  in  the  a-set  structural  coordinates, 
satisfy  the  eigenvalue  problem 

[Kaa]  [<i>ai\  =  [M^a]  [M  [0]  (2.1) 

where  [fl]  is  a  diagonal  matrix  of  the  rij  eigenvalues.  If  static  Guyan  reduction  is  applied, 
the  modes  can  be  recovered  to  the  f-set  by 

W  (2.2) 

Recovery  of  single-  and  multi-point  constrained  displacements  leads  to  the  g-  set  normal 
modes  [^^]. 

The  basic  assumption  of  the  modal  approach  is  that  the  structural  displacements  during 
structural  response  to  external  excitation  can  be  adequately  expressed  as  a  linear  combination 
of  the  baseline  modes: 

{Ua}  =  W{e}  (2.3) 

where  {uo}  is  the  a-set  structural  displacement  vector  and  {^}  is  the  vectors  of  generalized 
displacements. 

2.2  Modal  Mass  and  Stiffness  Matrices 

The  discrete-coordinate  g-set  stiffness  and  mass  matrices  are  assembled  in  the  design  loop 
by  adding  the  contributions  of  the  ridv  global  design  variables  to  the  contribution  of  the 
elements  that  are  not  subject  to  changes.  The  assembly  equations,  (5-43)  and  (5-44)  of  Ref. 
9,  when  all  the  structural  elements  are  linear  with  respect  to  the  design  variables,  can  be 
written  as 

[K,,]  =  [DKV]o  +  £  Vi[DKY]i  (2.4) 

»=i 
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where  [DKV]j  is  a  regular  linear  contributions  of  the  design  variable  Vj.  The  ma.ss  matrix  is 
assembles  by 

^dv 

[M,,]  =  [DMVjo  +  E  Vi[DMY]i  (2.5) 

t=i 

The  coefficient  matrices  in  Eqs.  (2.4)  and  (2.5)  are  stored  in  the  standard  ASTROS  data 
base  for  subsequent  sensitivity  analyses. 

The  premultiplication  of  Eq.  (2.4)  by  [<f)gi]'^  and  the  postmultiplication  by  yield  the 
assembled  generalized  stiffness  matrix 


^dv 

[Kii\  =  [DGKVjo  +  £ui[DGKV]i 

i=l 


(2.6) 


where 

[DGKVli  =  W’'[DKV]i[.^J  (2.7) 

The  rigid-body  related  portions  of  [iiiii]  are  zero.  The  premultiplication  of  Eq.  (2.5)  by  \4>g^ 
and  the  postmultiplication  by  yield  the  generalized  mass  matrix 


^dv 

[Mu]  =  [DGMV]o  +  £;«4DGMV]i 

»=1 


(2.8) 


where 


[DGMVJi  = 


(2.9) 


A  major  advantage  of  the  modal  approach  is  that  we  can  modify  the  generalized  stiffness 
and  mass  matrices,  due  to  design  changes,  and  repeat  the  analysis  without  returning  to  the 
large-size  finite-element  model.  When  all  the  design  variables  are  at  their  baseline  values 
both  [ATji]  and  [Mjj]  are  diagonal.  They  are  generally  full  in  the  modified  structure  (except 
for  the  rigid-body  rows  and  columns  in  [Ku]  which  remain  zero). 

Denoting  these  baseline  matrices  by  subscript  6,  we  can  rewrite  Eqs.  (2.6)  and  (2.8)  as 

[Ku\  =  [Ku\,  +  J2{vi  -  U6j[DGKV]i  (2.10) 

i=i 
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and 


^dv 

[Mii\  =  +  Y^ivi  -  ViJ[DGMV]i  (2.11) 

t=i 

Aeroservoelastic  analysis  with  moving  control  surfaces  requires  the  mass  coupling  [M^c] 
between  the  modal  coordinates  and  the  control  surface  deflection  modes.  For  the  baseline 
structure, 

Tldv 

=  [DGMVC]o  +  £  vb.  [DGMVCji  (2.12) 

i=l 

where 

(DGMVCIi  =  IDMVJifc]  (2.13) 


where  [4>g^  is  the  matrix  of  g-set  control  modes.  Each  control  mode  contains  the  kinematic 
structural  displacements  due  to  a  unit  deflection  of  the  associated  control  surface.  The 
modified  mass  coupling  matrices  is  and 

[Mi,]  =  [Mi,]b  +  £(ui  -  UbJpGMVCii  (2.14) 

i=l 


2.3  Normal  Modes  of  the  Modified  Structure 


The  Tti  normal  modes  [<^ai]  of  the  modified  structure  satisfy  the  eigenvalue  equation 

[KaaMai]  =  [Maa\[km  (2.15) 

where  [Q]  is  a  diagonal  matrix  of  the  corresponding  eigenvalues. 

It  is  assumed  that  the  normal  modes  of  the  modified  structure  are  linear  combinations 
of  the  baseline  modes  [4>]  at  all  the  discrete  coordinate  set  levels.  At  the  a-set  level  it  reads 

fc]  =  [4>aM  (2.16) 

where  [ip]  is  a  square  non-singular  matrix.  The  substitution  of  Eq.  (2.16)  in  Eq.  (2.15)  and 

premultiplication  by  yields  the  eigenvalue  problem 

[Kii][iP]  =  (2.17) 
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where  [Ki^  is  defined  in  Eq.  (2.6)  and  [M^i]  in  Eq.  (2.8).  Eq.  (2.17)  can  be  solved  for  the 
eigenvalues  [12]  and  the  eigenvectors  [•0] ,  with  [^]  normalized  such  that 

-  [/]  (2.18) 

which  yields 

m  (2.19) 

To  allow  the  formulation  below  to  be  applicable  to  the  baseline  structure  as  well  as  the 
modified  one,  we  define  the  baseline  [ip]  as 

The  application  of  Eq.  (2.17)  to  the  j-th  eigenvector,  {ipj},  its  differentiation  with  respect 
to  a  design  variable  Uj,  and  premultiplication  by  {ipj}^  yield  the  sensitivity  of  the  j-th 
eigenvalue, 

^  =  {’PiV  ([DGKVj.  -  n,[DGMV]j)  {^,}  (2.20) 

which  can  be  used  in  calculating  the  frequency  constraint  sensitivities.  In  order  to  apply 
Eq.  (2.20)  and  the  generalized  formulation  below  to  the  baseline  structure  as  well  as  to  the 
modified  ones,  the  baseline  modes  are  normalized  for  unit  generalized  masses  by  Eq.  (2.16) 
with 

(2.21) 

2.4  Generalized  Dynamic  Matrices 

2.4.1  Dynamic  equation  of  motion 

The  generalized  dynamic  matrices  are  constructed  in  the  h-set  coordinates  which  are  based 
on  the  i-set  modes  plus  extra  (e-set)  user-defined  coordinates.  As  mention  above,  the  formu¬ 
lation  below  assumes  that  there  are  no  extra  degrees  of  freedom,  such  that  =  nj.  However, 
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while  i-set  matrices  are  related  in  this  document  to  the  baseline  modes  [^ot])  h-set  matrices 
are  related  to  the  modified  modes 


[4>ah]  =  =  [<t>aM]  (2.22) 

where  [^ai]  is  kept  fixed  and  [^o^]  varies  during  the  optimization  according  to  [•0]  .  Another 

difference  between  the  h-set  and  the  i-set  matrices  is  that,  like  in  the  standard  ASTROS, 

the  h-set  matrices  include  the  effects  of  direct-input  structural  matrices  [Kgg]  and 

[Bgg]  which  are  not  included  in  the  eigenvalue  analysis  and  are  not  subject  to  changes  in  the 
design  process. 

The  time-domain  matrix  equation  of  motion  in  modal  coordinates  is 

[Mw.){|}  +  («**]{?}  +  [JTuKf}  =  {«,(()}  (2.23) 

where  the  generalized  mass  and  stiffiiess  matrices  are 

[ir«.)  =  (w  +  ra)  m  =  in]  +  W  (2.24) 

and 

[Mkkl  =  M’-  ((Miil  +  (M?))  M  =  [/]  +  m  (2.25) 

where 

The  generalized  damping  matrix  is  based  on  the  regular  damping  options  in  ASTROS, 
with  some  modifications, 

=  [9(u.»)l  W  +  [V'nBaW  +  .^[S]  (2.26) 

where 

m = 


8 


and  [a>/,]  =  The  optional  user-input  damping  parameters  here  are  those  used  in  the 

standard  ASTROS,  namely  the  general  structural  damping  g,  uj^  which  defines  the  equivalent 
viscous  damping,  and  the  modal  damping  table  g{ufh)- 

The  excitation  vector  in  Eq.  (2.22),  {Pfc},  can  contain  prescribed  external  inputs,  aero¬ 
dynamic  forces  which  are  discussed  in  the  following  section,  and  inertial  forces  due  to  control 
conmiands.  The  generalized  forces  due  to  prescribed  external  excitation  are 

{pn = iM^pm 

The  inertial  forces  due  to  control  surface  accelerations  are 

{Pl}  = 

where 

=  I'I’P  (IMJ  +  [JI41) 

where 

[Ml] = 

2.4.2  Sensitivity  of  Dynamic  Matrices 

The  change  of  modal  coordinates  (by  changing  [V’])  in  each  design  iteration  complicates  the 
sensitivity  analysis  compared  to  that  of  the  fixed-basis  analysis.  To  allow  simple  expressions, 
the  sensitivity  analysis  in  each  iteration  assumes  that  the  dynamic  matrices  of  the  next 
iteration  will  be  based  on  the  current  modal  coordinates.  Consequently,  the  differentiation 
of  [iiOih]  and  [Mhh]  with  respect  to  the  design  variables  are  based  on  Eqs.  (2.25)  and  (2.24) 
with  fixed  [•0],  and  [M?],  and  variable  {Kii\  [Mu].  The  resulting  sensitivities  are 

=  [l(-nDGKV|iW  (2.30) 

and 

^  [Mhh  Mfc,]  =  [0]^[DGMVi0  DGMVCi]  (2.31) 


(2.27) 

(2.28) 

(2.29) 
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2.4.3  Integrated  section  loads 

Section  loads  due  to  aeroelastic  response  are  calculated  by  the  summation-of-force  approach 
which  adds  the  effects  of  aerodynamic,  inertial  and  prescribed  forces.  The  effects  of  the 
prescribed  forces  are 

{PD  =  [uYipm  (2.32) 

where  [(figi]  is  the  matrix  of  integration  load  modes.  Each  load  mode  contains  the  kinematic 
structural  displacements  due  to  a  unit  displacement  of  a  reference  degree  of  freedom  defined 
by  the  user.  The  inertial  parts  of  the  section  loads  can  be  calculated  by 


{Pi} = 


(2.33) 


where  the  mass  matrices  are  first  defined  by 


Mia  Mi^  =  (.^51.)^  {\M„\  +  [JMj,])  [  ] 


(2.34) 


and  then  updated  by 

[Mu  Mlc]  =  [  Mu  Mu\  +  £(vi  -  U6i)  [  DGMVLIi  DGMVLCi  ]  (2.35) 

»=1 

from  which  is  calculated  by 


[Mlh]  =  [MuM 


(2.36) 


Following  the  sensitivity  formulation  above,  the  derivatives  of  the  mass  matrices  associated 
with  the  section  loads  are 


dvi 


Mlh  Mu]  =  [  DGMVLIi^  DGMVLCi  ] 


(2.37) 
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2.5  Aerodynamic  Forces 

The  unsteady  aerodynamic  routines  in  ASTROS  assume  that  the  structure  undergoes  har¬ 
monic  oscillations.  Complex  aerodynamic  force  coefficient  (AFC)  matrices  (at  a  given  Mach 
number)  for  n/  user-defined  ’’tabulated”  values  of  the  reduced  frequency  k  =  ufb/V  where  u> 
is  the  frequency  of  oscillations  and  6  is  a  reference  semichord.  These  complex  AFC  matrices, 
are  stored  in  a  data  base  for  subsequent  analyses.  Their  dimension  is  n*  x  n*. 
where  n*.  is  the  number  of  aerodynamic  degrees  of  freedom.  The  ZAERO  unsteady  aero¬ 
dynamic  module  also  calculates  the  coefficient  matrices  [Qkciiki)],  [QkG{i^i)]i  [QLk{iki)], 
and  [Qwi^ki)]  in  the  following  generalized  force  expressions.  These  matrices  are 
independent  of  the  structural  properties. 

The  generalized  unsteady  aerodynamic  forces  acting  on  the  structural  modes  of  a  linear 
aeroelastic  system  can  be  expressed  in  the  frequency  domain  as: 

{P^(ia;)}  =  -q[Qfcfc(tA:)]{^(ict;)}  -  g[<5hc(*fc)]{4(j‘«^)}  -  ^[QhG{ik)]{wa{w)}  (2.38) 

where  q  is  the  dynamic  pressure,  V  is  the  true  air  velocity,  and  {^}  ,  {(Jc}  and  {iog}  are 
the  vectors  of  generalized  structural  displacements,  ric  control  surface  commanded  deflec¬ 
tions  (actuator  outputs),  and  n©  gust  velocities.  [Qfch],  [<5fcc]  and  [Qho]  are  the  associated 
generalized  unsteady  aerodynamic  force  coefficient  matrices  calculated  by 

[  Qhh  Qhc  QhG  ]  =  [  Qkk<f>kh  Qkc  QkG  ]  (2.39) 

where 

[<l>kh]  =  [Gkf][(f>fh] 

where  [(?*/]  is  the  spline  matrix. 

The  unsteady  aerodynamic  section  loads  are 

{Pliiui)}  =  -g[QLfc(*fc)]{^(w)}  -  q[QLc{ik)]{Sc{ioj)}  -  ^[QLG{ik)]{wGiiu;)}  (2.40) 
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where  [Qx,)i]  is  calculated  by 


[Qlh]  =  [<?Lfe][^jb/i]  (2-41) 

Once  calculated  for  the  baseline  structure  with  the  original  modes  the  generalized 
force  coefficient  matrices  are  updated  during  the  optimization  by 

Qhh  Qhc  QhG 
Qlh  Qlc  QlG 


ip'^QiG 

Qia^  Qlc  Qlg 
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Chapter  3 

Rational  Aerodynamic 
Approximations 


Rational  aerodynamic  approximations  (RAA)  are  required  in  order  to  cast  the  dynamic 
aeroelastic  equation  of  motion  in  a  state-space  form.  The  approximation  options  described  in 
this  chapter  should  be  selected  for  best  combination  of  accuracy  and  efficiency  in  subsequent 
stability  and  response  analyses. 

The  approximations  are  constructed  for  the  baseline  model  with  the  aerodynamic  data 
associated  with  the  original  modes  Nevertheless,  the  subscripts  of  the  generalized 

aerodynamic  matrices  in  this  chapter  are  denoted  by  h  instead  of  i  for  consistency  with  the 
notations  in  other  documents.  The  resulting  approximation  matrices  remain  unchanged,  ex¬ 
cept  multiplication  by  normalization  matrices  (see  Chapter  4),  throughout  the  optimization 
process,  in  all  the  ASE  disciplines. 


3.1  Minimum-State  Approximation  Formula 


The  generalized  unsteady  aerodynamic  forces  in  Eqs.  (2.38)  and  (2.40)  can  be  expressed  in 
the  Laplace  s  domain  as: 


Qhh{s) 

QLh{s) 


Qhc{s) 

Qlc{s) 


QhG{») 

Qlg{s) 


(3-1) 
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The  AFC  matrices  are  normally  not  available  as  explicit  functions  of  s.  In  order  to  cast 
the  aeroelastic  equation  of  motion  in  a  time-domain  constant  coefficient  equation,  the  AFC 
matrices  have  to  be  described  as  rational  functions  of  s.  The  most  general  rational  function 
approximation  of  the  merged  AFC  matrices 


Nl  _  Qh  _  Qhh  Qhc  QhG 

_Ql  _  Qlh  Qlc  Qlg 

that  leads  to  a  state-space  aeroelastic  model  can  be  cast  in  the  form 

A/|0  ■A/,2  2  I  /rri  r  e>1\ — 1  r  nrl  /n 

where  p  is  the  nondimensional  complex  Laplace  variable  p  =  sb/V  and  all  the  matrix  coef¬ 

ficients  are  real  valued  [2].  The  common  [/2]  and  [^]  facilitates  the  formulation  of  section 
loads  as  output  variables  in  the  state-space  aeroelastic  equations  of  motion.  The  number 
of  structural  states  in  the  resulting  state-space  model  is  2n/^.  The  number  of  aerodynamic 
states  (n^)  is  equal  to  the  order  of  the  aerodynamic  root  matrix  [R]. 

Since  the  aerodynamic  data  is  given  for  harmonic  oscillations,  the  approximation  process 
starts  with  the  replacement  of  p  in  Eq.  (3.2)  by  ik^  where  k  is  the  nondimensional  frequency 
ujbfV .  Least-square  procedures  are  then  used  to  calculate  the  approximation  coefficients 
that  best  fit  the  tabulated  [Q{iki)]  matrices. 

Roger  [1]  dealt  with  [Q}^  only.  He  approximated  each  aerodynamic  term  [Qh{ik)]  sepa¬ 
rately,  but  with  ui  common  aerodyn2miic  roots. 


[<3k(ifc)]  =  [>1m1  +  iilA..]  -  +  Y. 

1=3  **  I  7i-2 

which  can  be  cast  in  the  form  of  Eq.  (3.2)  with 

[71-^  1  r  A«  ■ 

[£>a]  =  [/  /  ...],  [R]  =  -  72/  ,  [E]=  Am 

which  implies  that  the  resulting  number  of  aerodynamic  states  is  na  =  ni  x  n/,. 


(3.3) 


(3.4) 
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Elquation  (3.4)  shows  that  Roger  treated  the  elements  of  \E]  as  free  variables  in  the 
approximation  process,  while  those  of  [D/,]  were  fixed,  and  that  [i2]  has  repeated  roots.  The 
MS  method  [2,  3]  was  based  on  the  realization  that  the  number  of  aerodynamic  states  per 
desired  accuracy  can  be  reduced  significantly  by  treating  all  the  elements  of  both  [jDfc]  and 
[.E]  as  free  variables,  and  by  letting  the  diagonal  [jR]  be  with  distinct  negative  roots.  The 
number  of  MS  distinct  roots,  nj,  and  their  values  are  defined  by  the  analyst.  This  number  is 
typically  larger  than  that  in  Roger’s  approximation,  but  the  number  of  resulting  aerodynamic 
states,  Ua  =  ni  is  much  smaller.  Typical  MS  applications  [4-6]  used  4  to  6  roots  with  the 
span  of  their  absolute  values  similar  to  that  of  the  tabulated  reduced  frequencies.  The  MS 
applications  showed  small  sensitivity  to  the  root  values. 

3.2  Weighted  Least-Square  Fit 

To  facilitate  real- valued  algebra,  the  complex  approximation  expression  in  Eq.  (3.2),  with  p 
replaced  by  ik,  is  separated  into  real  and  imaginary  parts.  The  real  part  of  the  h  partition 
is 

|ft(fc)l  =  [Am]  -  P(Am]  +  (3.5) 

where 

lAr(fc)j  =  (4^/1  + 

is  diagonal,  and  the  imaginary  part  is 

[Gh(A:)]  =  k[AHi]  -  k[DH][K{k)][R][E]  (3.6) 

The  L  partition  of  Eq.  (3.2)  is  expressed  by  Eqs.  (3.5)  and  (3.6)  with  subscript  h  replaced 
by  L.  It  can  be  noticed  that  both  [Qh]  and  \Ql[  use  the  same  [R],  [isr(fc)]  and  \E\  matrices. 
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The  comparison  of  Eqs.  (3.5)  and  (3.6)  with  the  real  and  imaginary  parts,  [i^(A;^)]  and 
[G(A:/)],  of  the  tabulated  AFC  matrices  [Q(iA:/)]  provides  an  overdetermined  set  of  approxi¬ 
mate  equations.  The  task  is  to  find  the  free  approximation  coefficients  in  [Aq],  [Ai],  [A2],  [D] 
and  [E]  that  minimize,  under  some  constraints,  the  total  least-square  approximation  error 

~  ~  {3-7) 

V  • 

where  Wijt  is  the  weight  assigned  to  the  ij-th.  term  of  the  Ath  tabulated  AFC  matrix.  The 
weight  options  in  Eq.  (3.7)  are  discussed  in  Section  3.5. 

For  numerical  efficiency,  it  is  desired  to  decompose  the  overall  least-square  problem  into  a 
sequence  of  small-size  standard  least-square  problems,  each  based  on  the  general  approximate 
equation 

[PF*],  «  [W*Ub*}t  f<yri  =  l,nk  (3.8) 

where  {®*}  is  a  subset  of  unknown  coefficients  which  are  uncoupled  with  others,  {6*}/  is 
extracted  from  the  tabulated  data,  and  [A*]/  is  a  function  of  fc/,  the  aerodynamic  root  values 
and  the  constraints  associated  with  unknowns,  as  detailed  below  for  the  various  techniques. 
The  weight  matrix  is  a  diagonal  matrix  whose  terms  are  the  weights  Wijt  associated 

with  data  terms  in  {6*}^.  The  weighted  least-square  solution  for  {s*}  is  obtained  by  solving 
the  symmetric  problem 

[G]{x*}  =  {6}  (3.9) 

where 

[Cl = 

i 

{*>} = 

/ 

Various  approximation  cases  are  presented  in  Sections  3.3  and  3.4  by  defining  their  sequence 
of  least-square  solutions,  the  unknown  coefficients  {x*}  in  each  solution,  and  the  associated 
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data  vectors  {6*}/  and  approximation  matrices  [A*]i  of  Eq.  (3.8).  The  weight  are  discussed 
in  Section  3.5. 


3.3  The  Unconstrained  Problem 


The  approximation  coefficients  of  Eq.  (3.2)  are  determined  by  solving  a  sequence  of  weighted 
least-square  problems.  The  coefficients  of  are  determined  by  either  Roger’s  or  the  MS 
procedures.  The  coefficients  of  [Ql]  are  then  defined  with  [E]  of  the  [Qfc].  The  diagonal  [R] 
is  defined  by  the  user  in  all  the  cases. 

3.3.1  Roger’s  approximation 

The  structure  of  Eq.  (3.3)  indicates  that  the  coefficients  of  Roger’s  approximation  can  be 
found  for  each  aerodynamic  term  separately  by  performing  rih  x  {uh  +  nc  -f-  uq)  solutions  of 
Eq.  (3.9)  with 


= 


1  0  ^  A 


4^ -J’ 

(3.10) 

where  the  number  of  unknowns  in  each  solution  is  nj  +  3.  The  right-side  terms,  and  Ghij 
are  the  real  and  imaginary  parts  of  the  {i,j)  term  of  [(5fc(iA:/)].  The  resulting  {x*}  vectors 
and  the  user-input  lag  term  are  used  to  define  the  coefficient  matrices  in  Eqs.  (3.3)  and 
(3.4). 

3.3.2  Minimum-state  procedure 

With  both  [£)/,]  and  [E]  in  Eq.  (3.2)  being  unknown,  the  MS  problem  for  a  given  [R]  is 
nonlineax.  It  is  solved  iteratively  by  starting  with  an  initial  guess  of  [E;,]  in  which  at  least 
one  term  in  each  row  and  each  column  is  nonzero.  For  a  given  [E/,],  Eqs.  (3.5)  and  (3.6) 


f 

1  {*'}  =  1 


{»•},  = 


imply  that  the  unknown  MS  coefficient  matrices  [j4ho]>  [-^w]  and  [£?]  can  be  calculated 

by  performing  n/,  +  Wc  +  ng  column-by-column  solutions  of  Eq.(3.9)  with 


,  ,  ^  I  0  -kll  klD^K(h)  _  I  _  f  Fi^(h)  \ 

‘  [o  hi  0  -hDHK{k,)R\'  ^  (’  '  G»,(ife,)  / 


(3.11) 

where  the  j  indices  relate  to  the  j’-th  columns  of  the  respective  matrices.  The  number  of 
unknowns  in  each  regular  solution  is  -f  nj.  The  calculated  [£^]  is  then  used  to  update 
[ilfco],  [-4^:],  and  [D^]  by  performing  Uh  row-by-row  solutions  of  Eq.  (3.9)  with 


[^1^  = 


I  0  -kjl  kjE^K{ki) 

0  kj  0  -kiE^K{kt)R  ’ 


{x*}  = 


(3,12) 


where  the  i  indices  relate  to  the  i-th  rows  of  the  respective  matrices.  The  number  of  unknowns 
in  each  solution  of  Eq.  (3.12)  is  3(n^+nc+n(3r)+n/.  The  least-square  solutions  with  Eqs.  (3.11) 
and  (3.12)  form  a.  D  E  D  iteration  which  is  repeated  until  convergence  is  obtained  or 
until  a  specified  maximum  number  of  iterations  is  reached. 

The  iterative  nature  of  the  MS  procedure,  and  the  relatively  large  number  of  unknowns 
solved  for  simultaneously  in  each  iteration  of  the  regular  unconstrained  problem,  require  a 
considerably  larger  computation  time  than  that  of  Roger’s  method.  A  major  reduction  in  the 
MS  computation  efforts  was  obtained  in  Refs.  [2-7]  by  applying  3  approximation  constraints. 
The  efficient  solution  of  Eq.  (3.9)  presented  in  Ref.  [8]  allows  similar  computational  savings 
without  having  to  apply  constraints.  It  can  be  observed  that  the  structure  of  the  least-square 
equation  (3.9)  with  the  coefficients  and  variables  of  Eq.  (3.11)  is; 


(3.13) 


Cn 

0 

Ciz 

C'i4 

f 

0 

C22 

0 

^24 

0 

C33 

^^34 

CL 

^24 

rtT 

^34 

<^44  . 

1 
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where  [C22],  [C'ls]  and  [C'ss]  are  diagonal.  This  structure  facilitates  an  efficient  solution 
of  Eq.  (3.13)  through  a  sequential  reduction  of  the  problem  size.  The  first  row  of  Eq.  (3.13) 
yields: 

{A„,}  =  [Cn]-‘  ({6,}  -  (C.3KA2J  -  (3.14) 

The  substitution  of  Eq.  (3.14)  in  rows  2  to  4  of  Eq.  (3.13)  yields: 

■  ^22  0  C24  ]  f  Aw,  1  f  ^  1 

0  ^33  C34  ^  Ah2j  *  =  <  ^  '  (3.15) 

L  CS  Cl  C44  J  I  Ej  J  U.  J 

where  [(722]  and  [C33]  remain  diagonal.  The  next  step  is  the  extraction  of  {Aw^}  and  {Aw,} 
from  the  first  two  rows  of  Eq.  (3.15),  and  their  substitution  in  the  third  row.  The  result  is 

[cmy  =  {5}  (3.16) 

where 


[C]  =  C„-ClCi^C„-ClO^^C3t 


{6}  =  h-Clfi^h-ClC^b, 


Equation  (3.16)  is  of  order  nj  only,  with  symmetric  [C].  The  diagonality  of  [(7ii]  and 
[Cis]  in  Eq.  (3.14),  and  [(722]  and  [C33]  in  Eq.  (3.15)  implies  that  the  reduction  of  {Aw^}, 
{Aw,}  and  {Aw,}  does  not  involve  matrix  inversion  and  can  be  performed  term  by  term. 
The  reduction  of  the  i-th  row  in  {Aw,},  {Aw,}  and  {Aw,}  affects  only  the  i-th  row  and 
i-th  column  in  [^44]. 

The  construction  and  solution  of  Eq.  (3.16)  are  repeated  until  the  entire  [.E]  is  found 
for  the  last  calculated  [Dh\-  For  each  {Ej}  we  can  use  Eq.  (3.15)  for  {Aw^.}  and  {Aw,}, 
and  then  Eq.  (3.14)  for  {Aw,}-  This  recovery  is  required,  however,  only  at  the  end  of  the 
iterative  process,  unless  we  want  to  monitor  the  approximation  error  of  Eq.  (3.7)  during 
the  iterative  process.  The  same  reduction  and  solution  process  can  be  performed  for  [Aw], 
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[•Afci],  [>l/.2]  and  [£);,],  this  time  row  by  row,  with  the  coefficients  of  Eq.  (3.12)  using  the  last 
calculated  [E]. 

The  sequence  of  least-square  solutions  is  such  that,  theoretically,  the  total  error  should 
never  increase  when  the  process  progresses.  However,  numerical  difficulties  may  cause  it  to 
increase  when  high-order  approximations  are  applied  with  insufficient  amount  of  data  [6]. 
As  demonstrated  in  the  numerical  example  of  Ref.  [8],  the  reduction  process  reduces  the 
solution  time  of  a  typical  unconstrained  problem  by  more  than  80%. 

3.3.3  Integrated  load  coefficients 

The  coefficients  of  \Qi^  in  Eq.  (3.2)  are  calculated  by  a  direct  least-square  solution  with  the 
[E\  of  the  [Q/j]  approximation.  Equation  (3.9)  is  solved  nx,  times  (once. for  each  required 
integrated  load)  with  the  coefficients  and  variables  of  Eq.  (3.12),  with  subscripts  h  in  {®*} 
and  {6*}/  replaced  by  L. 

When  Roger’s  approximation  is  used,  the  resulting  [Dl]  is  an  nx,  x  [uh  *  ni)  matrix  and 
[Qi\  is  approximated  with  [i?]  and  [£]  of  Eq.  (3.4).  With  MS  approximation  the  size  of  [£>x,] 
is  nx  X  ni  and  [Qx]  is  approximated  with  [R]  and  [E]  of  [Qfc]. 

3.4  Application  of  Constraints 

Reference  [5]  presented  the  various  options  in  the  3-constraint  MS  procedure,  in  which  3 
constraints  have  to  be  assigned  to  each  aerod3miamic  term  even  if  they  are  not  desired.  Here, 
when  they  axe  not  required,  it  is  still  often  desired  to  obtain  exact  fits  at  specified  reduced 
frequencies  or  to  null  out  some  coefficients.  The  most  frequently  used  constraint  is  a  match 
of  the  steady-aerodynamics  data  (at  A:  =  0).  An  imaginary-part  data-match  constraint  at  a 
k  close  to  0  yields  dGijfdk  at  fc  =  0,  which  affects  the  quasi-steady  aerodynamic  damping, 
to  be  equal  to  that  of  the  tabulated  data.  Data-match  constraints  at  higher  k  values  are 
sometimes  desired  to  increase  the  accuracy  of  anticipated  flutter  mechanisms.  Even  with 
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the  reduction  technique  shown  above,  up  to  3  constraints  can  be  applied  in  a  way  that  saves 
some  computation  time.  On  the  other  hand,  it  always  increases  the  total  fit  errors,  and 
sometimes  causes  imwanted  wiggling  of  the  resulting  curve  fits  near  the  data-match  points 
[6]. 

The  constraint  formulation  is  given  below  for  the  approximation  of  [Qh])  but  is  applica¬ 
ble  to  the  terms  of  [Ql{  as  well.  The  optional  constraints  that,  when  applied,  reduce  the 
approximation  problem  size  are: 

1.  Steady  aerodynamics  match  enforced  in  Eq.  (3.5)  by  setting 


AhQij  = 


(3.17) 


2.  Either  imaginary-part  match  at  a  nonzero  k  =  kg  enforced  in  Eq.  (3.6)  by  setting 

^1.,  =  +  DH.K{k,)RBi  (3.18) 

or 

Ahu,  =  0  (3.19) 

3.  Either  real-part  match  at  a  nonzero  k  —  kf  enforced  in  Eq.  (3.5)  by  setting 

i  -  F^,{k,))  +  D^K{k,)B^  (3.20) 

or 

=  0  (3.21) 

Each  optional  constraint  in  Eqs.  (3.17)  thru  (3.21)  can  be  applied  for  a  different  sub¬ 
set  of  aerod3mamic  terms,  not  necessarily  for  all  {i,j)  terms  uniformly  as  done  in  previous 
developments  [4-6].  Different  aerodynamic  terms  can  be  assigned  with  different  match  fre¬ 
quencies,  namely  different  kg  and  kf  values  in  Eqs.  (3.18)  and  (3.20),  while  some  of  them 
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are  constrained  by  Eq.  (3.17)  and  the  others  are  not.  A  definition  of  the  default  constraint 
set  is  requested  in  the  ASE  module.  Other  constraint  sets  that  replace  the  default  set  can 
be  assigned  to  selected  columns  of  [Q/,]  and  [Ql]. 

The  constraints  are  applied  in  the  approximate  equation  level,  Eq.  (3.8),  by  eliminating 
terms  from  {a:*}  and  changing  [A*]/  and  {6*}/  accordingly.  This  is  done  before  the  construc¬ 
tion  of  Eq.  (3.9)..  The  process  is  described  below  with  the  unconstreiined  {a:*},  \A*\i  and 
{6*}/  of  Eq.  (3.11)  or  Eq.  (3.12)  being  the  starting  point. 

The  application  of  the  steady  aerodynamics  constraints  of  Eq.  (3.17)  eliminates  the  i-th 
term  in  {Aw),},  and  the  associated  column  in  [A*]/,  and  replaces  F}^-{kt)  in  {b*}  by 

=  PfH.ikt)  -  FH,,iO)  (3.22) 

The  application  of  the  imaginary-part  constraint  of  either  Eq.  (3.18)  or  Eq.  (3.19)  elim¬ 
inates  the  i-th  term  in  {Ai^},  and  the  associated  column  in  [A*]/.  The  match  constraint  of 
Eq.  (3.18)  also  replaces  the  i-th  row  of  the  block  [-kiDhK{kt)R]  in  [A*]/  by  [-kiDhiKg{kt)R] 
where 

[Kg{k,)]  =  [K{k,)]  -  [Ar(fc,)]  (3.23) 

and  the  Gh^-{ki)  term  in  {6*}  is  replaced  in  this  case  by 

Gfc,,  (3.24) 

'^9 

The  application  of  the  real-part  constraint  of  either  Eq.  (3.20)  or  Eq.  (3.21)  eliminates 
the  i-th  term  from  {Ah2j},  and  the  associated  column  from  [A*]/.  The  match  constraint 
of  Eq.  (3.20)  also  replaces  the  i-th  row  of  the  block  [k}DhK{ki)]  in  [A*]i  by  [k^Dh^Kf{k^)] 
where 

[Kf{h)]  =  [Ar(fc^)]  -  [K{kf)\  (3.25) 
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(3.26) 


When  Eq.  (3.20)  is  applied  together  with  Eq.  (3.17),  Fhi^{kt)  of  {6*}  is  replaced  by 

-  •P'<Ny(0)  -  I  {FH.,(h)  -  Fh.,{0)) 

Otherwise,  when  Eq.  (3.20)  is  applied  but  Eq.  (3.17)  is  not,  Fhij(ki)  becomes 

Fk.,  {kt)  =  -  !«»,(%)  (3.27) 

and  term  (i,i)  in  [A*]/  is  replaced  by 


lii 


(3.28) 


The  application  of  constraints  to  the  coefficients  of  Eq.  (3.12)  is  similar  to  that  shown 
above  for  Eq.  (3.11).  When  the  three  data-match  constraints  of  Eqs.  (3.17),  (3.18)  and  (3.20) 
are  applied  simultaneously  to  all  the  (i,  j)  terms,  and  when  all  the  terms  are  assigned  with 
the  same  match  frequencies  kg  and  kf,  the  least-square  matrices  in  Eq.  (3.11)  are  reduced  to 


kjDhKfiki) 
-ktDhKg{ki)R  ’ 


(3.29) 


where  Fh^-{ki)  and  Ghij{ki)  axe  defined  in  Eqs.  (3.26)  and  (3.24)  respectively,  and  Eq.  (3.12) 
is  reduced  to 


k]E'^Kf{ki) 
-kiE^Kg{ki)R  ’ 


K}  =  {Dl}, 


(3.30) 


Since  constraints  are  applied  in  the  approximate  equation  level,  Eq.  (3.8),  and  the  re¬ 
duction  process  in  Eqs.  (3.13-3.16)  is  performed  in  the  solution  level,  all  constraints  have  to 
be  applied  before  the  reduction  starts.  It  is  clear  from  the  previous  section  that  constraints 
can  be  applied  selectively  with  each  aerodynamic  term  having  a  different  number  and  type 
of  constraints.  To  summarize  the  solution  process,  a  single  D  E  —¥  D  iteration  in  MS 
approximation  of  [Q^]  is  performed  in  the  following  steps; 
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1.  Construction  of  [j4*]j  and  {6*)-/  of  Eq.  (3.11),  with  selected  constraints  from  Eqs.  (3.17- 
3.21),  for  the  j  =  l  column  of  [Qh(p)]  in  Eq.  (3.2). 

2.  Construction  of  \C\  of  Eq.  (3.9)  for  this  column,  whose  terms  correspond  to  {Ejy  and 
to  the  terms  in  {^4^)^  },  {-4/,!,}  and  }  that  were  not  constrained  in  step  1. 

3.  Reduction  of  [C]  to  obtain  [C]  of  Eq.  (3.16),  and  solution  of  Eq.  (3.16)  for  {Ej}. 

4.  Repetition  of  steps  1-3  for  j  =  2,nh  +  ric  +  ng  to  calculate  the  entire  [E\. 

5.  Application  of  steps  1-4  to  the  calculation  of  row  by  row,  starting  with  the  coef¬ 
ficients  of  Eq.  (3.12). 

6.  Recovery  of  the  constrained  terms  in  [A;,o],  [A^i]  and  [A/,2]  by  Eqs.  (3.17-3.21),  and 
the  unconstrained  terms  by  Eqs.  (3.14)  and  (3.15). 

With  Roger’s  approximation,  steps  1-4  and  6  are  performed  once,  starting  with  Eq.  (3.10). 
To  approximate  [Qi]^  only  a  single  E  D  iteration  is  performed. 

3.5  Data  Weighting 

Least-square  curve  fitting  tends  to  produce  smaller  percentage  errors  at  data  points  of  large 
numerical  values.  However,  large  numerical  values  do  not  always  reflect  the  actual  impor¬ 
tance  of  accmate  fit  of  an  element  in  the  aerodyneimic  matrices.  To  eliminate  the  effect  of 
the  way  the  structural  normal  modes  are  normalized,  they  are  normalized  to  unit  gener¬ 
alized  mass  before  the  generalized  aerodynamic  matrices  are  calculated.  The  weighting  in 
the  least-square  problem,  Eq.  (3.8),  can  be  either  uniform,  with  all  Wiji  =  1,  or  based  on 
the  automatic  physical  weighting  process  [4]  whose  expressions  are  given  below.  The  phys¬ 
ical  weighting  can  be  applied  only  to  the  aerodynamic  terms  which  appear  in  the  dynamic 
equation  of  motion,  namely  [Qh]- 
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3.5.1  Measures  of  aeroelastic  importance 

The  physical-weighting  algorithm  developed  in  [4-6]  was  designed  to  weight  each  term  of 
the  tabulated  data  such  that  the  magnitude  of  the  weighted  term,  Eq.  (3.36),  indicates 
its  “aeroelastic  importance”.  The  idea  is  that  the  weight  assigned  to  a  data  term  should 
be  proportional  to  the  estimated  effect  of  a  unit  approximation  error  on  a  representative 
aeroelastic  property.  Different  representative  properties  are  selected  below  for  the  structural, 
control  and  gust-related  partitions  of  the  [Qhiiki)]  matrices.  The  error  effects  are  estimated 
by  the  differentiation  of  the  selected  aeroelastic  properties  with  respect  to  the  aerodynamic 
terms. 

The  measures  of  aeroelastic  importance  are  based  on  the  frequency-domain  equilibrium 
equation 

[C«.(ii)l{«<fc)}  =  -  9l0i«(i4)l)  Wi*!)}  -  (3.31) 

where 

1.2y2  ikVa 

[C'w.(ifc)]  =  +  [Kii\  -H  q[Qhh{ik)] 

The  determinant  of  [C'wi(ifc)]  is  sometimes  called  the  flutter  determinant.  It  becomes  zero 
when  the  dynamic  pressure  q,  the  true  airspeed  V  and  the  reduced  frequency  k  axe  at  their 
flutter  boundary  values.  The  values  of  V  and  q  reflect  design  flight  conditions  at  which  the 
open-loop  system  is  stable  for  all  the  tabulated  kt  values.  The  modal  damping  g  is  one  of 
the  user-input  weighting  parameters  discussed  in  Section  3.5.2. 

The  weights  assigned  to  the  terms  of  a  data  matrix  are  based  on  their  eff^ect 

on  the  determinant  of  the  system  matrix  of  Eq.  (3.31).  The  absolute  values  of  the 

partial  derivatives  of  this  determinant  with  respect  to  Qhhij{ikt),  divided  by  the  determinant 
itself,  is  shown  in  [4]  to  be  the  ij-th  term  of  the  weight  matrix 
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The  weights  assigned  to  the  terms  in  the  j-th  column  of  a  data  matrix  is  based 

on  the  open-loop  Nyquist  frequency  response  of  the  7-th  actuator  to  excitation  by  the  j-th 
control  surface.  The  actuator  response  is  related  to  the  modal  response  by 


(3.33) 


where  [T{iu^)]  is  a  matrix  of  transfer  functions  relating  actuator  outputs  to  sensor  inputs, 
and  is  the  matrix  of  modal  deflections  at  sensor  inputs*  The  physical- weighting  transfer 
functions  in  [T(ia;)]  should  be  simple  and  unrelated  to  particular  aeroservoelastic  parame¬ 
ters,  Structural,  narrow-band  filters  with  high  sensitivity  to  parametric  changes  should  not 
be  included  as  it  may  result  in  the  assignment  of  low  weights  to  important  aerodynamic 
data  terms.  Hence,  the  physical  weighting  algorithm  assumes  transfer  functions  which  are 
based  on  a  third-order  actuator,  multiplied  by  a  control  gain.  The  magnitude  of  the  partial 
derivative  of  the  Nyquist  signal  with  respect  to  Qhci^{ikt)  is  shown  in  [4]  to  be  the  ij-th  term 
of  the  weight  matrix 

=  q  ^  (3.34) 

The  weights  assigned  to  the  terms  of  the  j-th  column  of  a  data  matrix  [Qhaiih)]  is 
based  on  the  power  spectral  density  (PSD)  of  the  open-loop  response  of  a  selected  structural 
acceleration  to  continuous  gust,  derived  from  Eq.  (3.31), 


1  [^hfc(**/)]  ^{QhGj  {ikl)} 


(3.35) 


where  [^2/1^]  is  a  row  vector  of  modal  displacements  at  the  selected  response  point,  and 
is  the  PSD  function  of  the  associated  gust  velocity, 

^2  r  ^  .  n/l  r  /.\2 


TT  V[l  +  iktL,m^ 
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where  j-th  gust  RMS  velocity,  and  Lg  is  the  gust  characteristic  length.  The  partial 
derivative  of  yj (ki)  with  respect  to  QhOij {ikt)  is  the  ij-ih.  term  of  the  weight  matrix 

[Whg\  =  ^ \[4>.h][Chh{ikt)]-^f  (3.36) 

where  in  an  no  x  ng  diagonal  matrix  whose  elements  are 

3.5.2  Weight  scaling  and  modifications 

The  variations  of  terms  in  the  weight  groups  Eq.  (3.32),  [Whc]t,  Eq.  (3.34),  and 

[Wholh  Eq.  (3.36),  with  k  may  have  very  sharp  peaks.  In  addition,  the  peak  values  of  many 
terms  may  be  several  orders-of-magnitude  smaller  than  other  peaks.  The  extreme  weight 
variations  have  the  effect  of  neglecting  much  of  the  data,  which  may  cause  numerical  ill- 
conditioning  problems  and  unrealistic  curve  fits.  To  ensure  realistic  interpolation  between 
the  tabulated  k  values,  and  to  facilitate  the  application  of  the  resulting  aeroelastic  model 
to  a  variety  of  flow  conditions,  structural  modifications,  and  control  parameters,  it  may  be 
desirable  to  moderate  the  weight  variations.  This  can  be  done  by  one  or  a  combination  of 
the  following  means: 

•  assignment  of  a  relatively  large  damping  parameter  g  in  Eq.  (3.31), 

•  widening  the  weight  peaks, 

•  scaling  up  the  extremely  low  weights. 

Peak  widening  is  performed  in  cycles  applied  sequentially  to  each  weights,  W»j(A:/),  of 
the  3  weight  groups.  In  each  cycle,  Wij{kt)  is  changed  to  max{l^ij(A:/_i),  Wij{ki),  Wij{ki^i)} 
of  the  previous  cycle.  The  weight  matrices  are  then  normalized  and  combined  to  the  final 
weight  matrix 

[WH]t=[[WHH]t  [Whc]i  [WhG]i]  (3.37) 
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where  a  term  in  [Whh]/  is 


where 


WhH, 


'ijt 


=  ^max| 


W. 


cut 


m&XijiWhhi^}  Whh,^ 


Whku 


(3.38) 


Whhii  =  {|<3wk,  J  (3.39) 

and  the  terms  of  [W)ic]/  and  [PV/,g]/  are  calculated  similarly  (but  separately).  The  upscale 
parameter  Wcut  is  defined  by  the  user.  The  resulting  magnitudes  of  the  weighted  terms, 

=  (3.40) 

fall  between  Wcut  a^d  1.0  when  the  value  of  1.0  typically  appears  only  once  in  each  group.  The 
modified  physical  weighting  is  actually  a  compromise  between  the  the  unmodified  one  (with 
n^id  and  Wcut  equal  zero)  and  the  data-normcdization  weighting  of  Ref.  [4].  With  n^d  —  nt 
and  Wcut  —  1-0,  all  the  physical-weighting  effects  are  suppressed  and  the  weighting  becomes 
a  data-normalization  one.  Recommended  parameters  in  typical  cases  [5]  are  riyjd  =  2  and 
Wcut  =  0.01.  Vcirious  applications  demonstrated  that  the  resulting  aeroservoelastic  models 
were  adequate  for  analyses  with  large  variations  of  dynamic  pressures  [4-6,  11],  control  gains 
[12],  and  structural  parameters  [13]. 

3.5.3  Relative  importance  of  modes 

The  physical  weights  can  be  used  to  rate  the  vibration  modes  according  to  their  relative 
aeroelastic  importance.  Based  on  the  magnitudes  of  the  weighted  aerodynamic  data  terms, 
Eq.  (3.40),  calculated  with  n^d  and  Wcut  equal  zero,  three  modal  measures  of  aeroelastic 
importance  are  defined  for  each  structural  vibration  mode  by 

Qhi  =  max{|0h/K,(i*:r)|},  (?*  =  inM{|gfcc.,(ifc/)|}  , 

Qgi  =  max{l0fcGo(jfc/)|}  (3.41) 
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These  measures  can  be  interpreted  as  indicators  of  the  aeroelastic  activity  of  the  i-th  vibra¬ 
tion  mode,  on  a  scale  of  0  to  1,  in  three  categories:  a)  influence  on  the  open-loop  system 
roots  (Q*);  b)  role  in  the  a«roservoelastic  loop  (Q*);  and  c)  contribution  to  gust  response 
((5p.  Being  based  on  a  limited  analysis,  these  measures  should  be  used  with  caution.  Their 
main  usage  is  in  supplying  physical  insight  and  in  pointing  out  the  structural  modes  that 
can  be  eliminated  in  subsequent  analyses  from  the  model  without  causing  significant  errors. 
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Chapter  4 

State-Space  Aeroservoelastic 
Equations 


The  time-domain  ASE  model  for  stability  analysis  is  constructed  from  the  separate  models  of 
the  aeroelastic  plant,  the  sensing  and  actuation  models  and  the  control  system,  all  expressed 
in  state-space.  In  the  aeroelastic  model,  each  modal  coordinate  is  represented  by  two  states: 
the  modal  displacement  and  its  velocity.  Rational  approximation  of  the  unsteady  AFC 
matrices  in  the  Laplace  s  domain  facilitates  the  augmentation  to  incorporate  the  aerodynamic 
states  in  the  model.  The  system  matrix  of  the  plant  model  can  be  used  for  flutter  analysis  of 
the  open- loop  aeroelastic  plant.  The  control  system  includes  the  control  surfaces  driven  by 
actuators,  sensors  related  to  the  structural  degrees  of  freedom,  and  a  linear  MIMO  control 
law  that  relates  the  actuator  inputs  to  the  sensor  readings.  The  full  ASE  model  can  be  used 
for  control  margin  computation  and  for  closed-loop  flutter  analyses.  Since  only  stability  and 
flutter  issues  are  addressed  in  this  chapter,  no  external  inputs,  such  as  pilot  commands  and 
wind  gust  inputs,  are  incorporated  in  the  model. 

4.1  Aeroelastic  Model 

The  second-order  time-domain  equation  of  motion  of  the  structure  was  defined  in  Eq.  (2.23), 
with  the  generalized  external  forces  represented  by  {Pfc(t)}.  For  stability  analysis  of  the 
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open-loop  aeroelastic  system  or  the  closed-loop  aeroservoelastic  system,  {P/i(t)}  includes 
aerodynamic  forces  due  to  structural  dynamics,  and  inertial  and  aerodynamic  forces  due  to 
control-surface  motion. 

The  control-induced  inertial  forces  are  defined  in  Eq.  (2.28).  The  generalized  aerody¬ 
namic  forces  are  defined  in  the  frequency  domain  in  Eq.  (2.38).  The  aerodynamic  force 
coefficient  (AFC)  matrices  are  first  calculated  at  several  user-defined  tabulated  reduced  fre¬ 
quency  values,  fc/,  similar  to  a  regular  frequency  domain  flutter  analysis.  The  tabulated 
matrices  are  used  for  approximating  the  AFC  matrix  as  a  rational  function  of  k  in  the  entire 
frequency  domain,  as  described  in  Chapter  3.  An  expansion  to  the  entire  Laplace  domain  is 
performed  by  replacing  ik  in  the  rational  expression  by  the  non-dimensional  Laplaoe  variable 
p  —  shiv,  which  yields  Eq.  (3.2).  The  substitution  of  p  =  shfV  in  the  expression  for  [<5/i(p)] 
of  Eq.  (3.2)  gives 

[l5/.(«)]  =  l^wl  +  M  (w^  -  jM)  [£]»  (4.1) 

The  [A/u]  and  [E]  matrices  are  column  partitioned  as 

=  A^]  (n  =  0,1,2),  [B)  =  l«h 

As  stated  at  the  beginning  of  Chapter  3,  the  coefficient  matrices  in  Eq.  (4.1)  are  calculated 
for  the  baseline  structure  with  the  generalized  matrices  associated  with  the  original  modes 
before  the  transformation  of  Eq.  (2.21)  is  made.  The  actual  modes  that  serve  as 
generalized  coordinates  in  each  optimization  iteration  are  [0afe]  which  relates  to  via  [V’] 
using  Eq.  (2.22).  Consequently,  the  coefficient  matrices  of  Eq.  (4.1)  have  to  be  updated  in 
each  iteration  by 

[Ahhn  Ahc„]  =  [i>]'^[Aiini’  Aicn],  [Dh]  =  [^]^[A],  [Eh]  =  [Ei][if]  (4.2) 

where  the  subscripts  i  relate  to  the  matrices  calculated  in  the  rational  approximation  process, 
denoted  in  Chapter  3  with  subscripts  h. 


31 


To  facilitate  state-space  formulation,  an  augmenting  aerodynamic  state  vector  of  dimen¬ 
sion  Ua  is  defined  by  its  Laplace  transform  as 

=  ([I]s  -  |[/e])  ((Bki{f(*)}  +  [£;,i{t(»)}) »  (4.3) 

Equation  (3.1)  without  the  gust  term  and  Eqs.  (4. 1-4. 3)  yield  the  s-plane  generalized  aero¬ 
dynamic  force  vector 

{PM}  =  -  9[W^)]{e(5)}  -  9[Q;.c(^)]{<Jc(^)} 

=  —q  -I-  -I-  {^(5)} 

-  <1  ([^/4C„]  +  -f  ^[44,,]^^)  {4(s)} 

-  g[£>fc]{sa(s)}  (4.4) 


Equations  (2.23),  (2.28)  and  the  time-domain  version  of  Eq.  (4.4)  yield  the  state-space 
open-loop  aeroelastic  equation  of  motion 


{®oe}  —  [44ae]{2'ae}’  d"  [.®ae]{^oe}’ 


(4.5) 


where 


{aJae}  =  ^ 

1  ^  1 

^  {Uae}  =  •< 

ft) 

i  J 

1  ic  1 

[^e]  = 

[Bae]  = 


0  [/]  0 
-[M]-'  -[M]-'  [Bhh  +  ^44,^,]  -q[M]-\DH] 

0  f[i2] 

0  0  O' 

-q[M]-' [A;,,,]  -[M]-^[Mhc  +  ^Ahc,] 

0  [E,]  0 
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[m]  =  [Mhh]  +  ^ 

The  number  of  states  in  Eq.  (4.5)  is  2nh  +  ria- 

The  outputs  of  the  aeroelastic  plant  are  sensor  readings.  It  is  assumed  here  that  the 
sensors  measure  either  structural  displacements,  velocities  or  accelerations,  and  that  the 
measurements  axe  perfect.  Sensor  dynamics  can  be  modeled  in  series  with  the  sensors’ 
outputs  and  incorporated  in  the  ASE  model  as  part  of  the  control  system,  discussed  in  the 
sequel.  The  outputs  are  assumed  to  be  linear  combinations  of  the  structural  state  response. 
The  combinations  are  defined  by  the  modal  displacement  (or  rotation)  row  vector  [<f)y]  at  the 
sensor  location.  Thus,  a  displacement  sensor  reads 


yd=[<f>y  0  0 

1 

(4.6) 

A  velocity  (rate)  sensor  reads 

O 

O 

II 

{®ae} 

(4.7) 

and  an  accelerometer  reads 

Va  = 

=  -[4>y][M]-^{[[KHh  +  qAHho]  qlDn]  ]  {Xae}  + 

+  [  ^[-^fcco]  ]  {^ae})  (4.8) 

In  general,  the  output  of  n,  sensor  readings  can  be  expressed  by 

{Vae}  =  [C'ae]{a;ae}  +  [I>ae]{«ae}  (4.9) 

4.2  Actuator  Model 


The  dynamic  model  of  the  actuator  driving  the  i-th  control  surface  is  specified  by  a  transfer 
function  having  the  form 


<^ci(^)  ^  _ OiS _ 

Uaci  (s)  s^  +  anS^  +  0,23  +  0,3 


(4.10) 
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where  Uaa  is  the  servo-commanded  (actuator  input)  control  surface  deflection.  The  DC 
(5  =  0)  gain  of  the  transfer  function  in  Eq.  (4.10)  is  one,  and  thus  in  steady  state  Uaa  =  S^. 

The  order  difference  of  three  between  the  numerator  and  denominator  in  Eq.  (4.10)  is 
justified  on  physical  groimds:  second  order  difference  represents  the  control  surface  deflec¬ 
tion  response  to  a  force  input  and  a  first  order  lag  used  to  model  the  realistically  limited 
bandwidth  of  the  actuator.  Due  to  this  order  difference,  6i,  Si  and  Si  can  be  defined  as 
independent  states  in  the  actuator  state-space  model,  and  thus  used  directly  as  inputs  to 
Eq.  (4.5)  and  in  connections  to  acceleration  sensors.  Higher  order  actuator  dynamics  can  be 
defined  by  connecting  additional  transfer  functions  in  series  to  Uaa  ■  These  additional  trans¬ 
fer  functions  can  be  included  as  part  of  the  control  system  model  discussed  in  the  following 
section. 

A  state-space  realization  of  the  actuator  dynamics  of  Eq.  (4.10)  is 


where 


■  0 

1 

0  ■ 

^  1 

Xaci}  = 

0 

0 

1 

{®ac.}  +  ^ 

0 

— Oi3 

— ®»2 

—Oil  _ 

1  J 

(4.11) 


< 


4 


1 


For  system  with  nc  >  1  actuators,  the  state- space  model  of  all  the  actuators  is  arranged 
so  that  the  total  actuator  state  vector  {xac}  equals  to  the  input  vector  {u^e}  of  Eq.  (4.5). 
Thus  the  actuator  state  vector  and  inputs  are 


f  f  1 

{Sac}  =  < 

r 

lij 

II 

: 

where  {5^}  ==  [^d  ^C2  *  *  *  The  state-space  equation  for  {xac}  is 


{®ac} 


fl2neXnc 

A. 


l2ncX2nc 


{Sac}  + 


Osnc  Xtic 

A. 


^OCJ 


Kc} 


(4.12) 
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where  Oanexne  l2ncX2ttc  respectively,  zero  and  identity  matrices  of  appropriate  dimen¬ 
sions,  and  [Aacj] ,  *  =  1,2, 3  are  diagonal  matrices  defined  by 

[i4acj]  —  diag{flij,  •  •  • ,  Ojic*} 

Equation  (4.12)  can  be  stated  in  a  compact  form  as 

{Xac}  =  Kc]  {lac}  +  [Bac]  M  (4.13) 

Since  {ajac}  =  {«ae},  the  augmentation  of  the  aeroelastic  states  {sae}  of  Eq.  (4.5)  to 
include  the  actuator  states  {sac}  of  Eq.  (4.13)  yields 

{ip}  =  [Ap]{xj,}  -h  [Bp\{up}  (4.14) 

where 


The  output  (sensor  measurement)  equation  (4.9)  becomes 

{Vp}  =  [<^p]{®p}  (4-15) 

where 

[C7p]=[aae  Doc] 

Note  that  in  equation  (4.15)  there  is  no  direct  feed-through  between  the  input  {up}  and 
the  output  {j/p}.  This  results  from  the  limited  bandwidth  of  the  actuators. 

4.3  Control  System  Model 

The  control  system  is  modeled  as  an  interconnection  of  three  types  of  basic  control  elements 
in  addition  to  a  variable  gain  matrix.  The  interconnections  within  the  control  elements  and 
between  them  and  the  aeroelastic  (AE)  system  can  be  either  fixed  or  through  the  variable 
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Control  System 

Figure  4.1:  The  control  system  interconnection  model.  Thick  lines  represent 
fixed  connections,  yp^^  i  =  1,2,3  are  the  sensor  reading  inputs 
and  Up-,i  =  1, . . . , 5  are  outputs  to  the  actuators. 

gain  matrix,  the  elements  of  which  are  subject  to  subsequent  stability  margin  analysis,  as 
demonstrated  schematically  in  Figure  4.1. 

The  ASE  model  can  be  used  for  parameter  optimization  of  the  AE  system.  In  this 
optimization  process,  the  control  system  is  assumed  to  be  constant,  except  of  possible  vari¬ 
ations  in  the  variable  gain  matrix  elements.  Thus,  only  the  AE  and  variable  gain  matrix 
are  changing  in  the  ASE  model  during  the  optimization.  To  avoid  unnecessary  recalculation 
of  the  entire  ASE  model  due  to  those  changes,  the  control  system  with  the  fixed  connec¬ 
tions  is  constructed  only  once  and  is  reused  in  the  ASE  model  reconstruction  throughout 
the  optimization  process,  thus  reducing  computational  load.  Clearly,  if  the  control  system 
is  updated  as  a  result  of  significant  changes  in  the  AE  model,  the  entire  ASE  model  has  to 
be  reconstructed. 
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4.3.1  Control  System  Elements 

The  control  system  is  constructed  from  three  types  of  elements:  single-input-single-output 
(SISO)  elements  specified  as  transfer  functions,  multi-input-multi-output  (MIMO)  elements 
defined  in  state-space,  and  zero-order  junctions.  The  parameters  of  these  elements  are  as¬ 
sumed  to  be  constant.  The  overall  control  system  variations  are  accounted  in  the  above 
mentioned  variable  gain  matrix. 

SISO  Control  Elements 


A  SISO  control  element  is  defined  by  a  proper  transfer  function  Tae(s),  which  is  then  realized 
in  state-space  inside  the  ASE  module.  A  general  transfer  function  from  the  input  Wae,»  to 
the  output  y,e^i  is  expressed  by  a  ratio  of  polynomials  in  s 


rp  /  \  _  yse,i{s)  _  -f  hx^jS^  ^  +  •  •  •  +  &n,i 

- -H  an,i 


The  controller  canonical  realization  of  Eq.  (4.16)  is 


yse,i  ~ 


where 


[Cae.i] 


f(n-l)x(n-l) 


{Bge^i} 


=  6o,t 


(4.16) 


(4.17) 


(4.18) 


The  number  of  states  in  {ause,*}  is  equal  to  the  order  n  of  the  denominator  polynomial  in 
Eq.  (4.16). 
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MIMO  State-Space  Control  Elements 


A  MIMO  state-space  control  element  is  defined  by  its  order  (number  of  states),  the  number 
of  inputs  and  outputs,  and  the  associated  state  space  dynamics  and  output  matrices.  These 
matrices  are  defined  by  the  user  as  direct  matrix  input  (DMI)  entries  or  by  INPUTT2 
MAPOL  commands.  This  facilitates  interaction  with  external  control  design  codes.  The 
equations  of  a  MIMO  state-space  element  are 

{®Tne,»}  =  [-^e,i]{®me,»}  +  [■5me,t]{^me,»} 

(4.19) 

A  special  case  of  Eq.  (4.19)  is  a  SISO  element  specified  in  state-space  format,  in  which  case 
[5TOe.i]  is  a  column  vector,  [Cme,.]  is  a  row  vector  and  is  a  scalar. 

Zero-Order  Control  Elements 

A  zero-order  element  is  defined  by  its  number  of  inputs  and  number  of  outputs.  This  is 
a  junction  element  in  which  every  output  is  a  weighted  sum  of  the  inputs.  The  element 
equation  is  simply 

{j/je,*}  =  [■Ojei»]{^ie,»}  (4.20) 

where  the  element  of  are  the  various  weights,  which  can  be  positive,  negative  or  zero. 

Clearly,  the  zero-order  element  is  a  special  case  of  a  MIMO  element.  It  is  separately 
introduced,  however,  to  simplify  the  input  format  and  internal  assembling. 

4.3.2  Fixed  Connections  in  the  Control  System  Model 

The  overall  control  system  is  constructed  as  an  interconnection  of  the  above  defined  control 
elements  and  the  variable  gain  matrix,  discussed  in  more  detail  in  the  next  section.  To 
construct  this  interconnection,  all  the  control  system  elements  are  combined  into  one  (not 


38 


connected  yet)  state-space  model 


<  Vme  > 

I  Vie  J 


A,,  0 

0 

Cse  0 
0  Gme 
0  0 


^se 

^me 

Xme 


+ 

Bse 

0 

L 

0 

Bme 

o  o 

1 _ 1 

1 

Use  1 

■  Dse 

0 

0  ■ 

f  Use  1 

+ 

0 

Dme 

0 

< 

Ume 

0 

0 

1 

[  J 

(4.21) 


where  the  various  blocks  represent  all  the  elements  of  the  same  type.  For  each  block,  its 
state,  input  and  output  vectors  combine  all  the  element  vectors,  and  the  state-space  matrices 
are  block  diagonal.  As  an  example,  the  SISO  elements  are  grouped  as  follows: 


Use,l 

\  S/»e.l  \ 

{®*e}  =  ■< 

^  ^se.nse  J 

►  {tise}  = 

^  Use,n,e  , 

^  {y>e}  =  < 

: 

[  j 

—  diag  [  i4ge,l  ■  ■  ■  Ase.tije  ]  I-®**]  —  diag  [  ■  ■  ■  Bte,n,e  ] 

[C„]  =  diag  [  Cse,l  •  •  •  C'se.n.e  ]  l^.e]  =  diag  [  Dsc,l  ■  ■  ■  ] 

where  rise  is  the  number  of  SISO  elements.  The  same  construction  is  performed  for  the  71^6 

MIMO  and  the  Uje  zero-order  elements. 

For  simplicity,  it  is  assumed  that  the  input  and  output  of  a  fixed  connections  cannot 
appear  in  another  connection.  Thus,  if  an  output  of  an  element  is  used  as  an  input  to  more 
than  one  other  element,  it  should  be  first  connected  to  a  “splitting”  zero-order  junction  with 
unity  weights,  the  outputs  of  which  can  then  be  used  as  necessary.  Similarly,  if  an  input 
of  an  element  is  a  sum  of  outputs  of  several  other  elements,  it  should  be  preceded  with  a 
summing  junction  (many  inputs  and  one  output).  To  avoid  confusion,  an  output  of  a  control 
element  cannot  be  connected  to  its  own  input.  After  the  fixed  connections  are  applied,  the 
associated  inputs  and  output  are  eliminated  from  the  equations  of  motion. 

In  a  fixed  connection  among  control  elements,  an  input  of  a  control  element  Ug-  is  set 
equal  to  an  output  of  another  control  element  j/g,  ,  namely  —  j/e,  *  To  perform  all  these 
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connections,  we  rewrite  the  state-space  model  of  Eq.  (4.21)  and  re-arrange  its  inputs  and 
outputs  as  follows: 


Xe  ^  A-e  j  Xg  ^  -|- 


Vet 

y^2 


a: 

a. 


D, 

D, 


eii 


*21 


(4.22) 


where  {uej}  and  {yei)  are  the  inputs  and  outputs  which  are  involved  in  fixed  connections. 
Obviously,  the  dimensions  of  {tiej}  and  {yej}  are  equal.  The  assumption  that  and  ye-  of 
a  certain  fixed  connection  cannot  appear  in  another  connection  implies  that 


KJ  =  [/2]{yeJ 


(4.23) 


where  [/2]  is  square  matrix  with  ones  at  the  (i,  j)  connection  entries  and  zero  elsewhere,  i.e., 
each  row  or  column  of  [I^  has  only  one  non-zero  entry  set  to  one.  Substituting  Eq.  (4.23) 
in  Eq.  (4.22)  yields  the  control  system  equation 


{Xc}  —  [>lc]{®c}  +  [•®c]{tic} 

{yc}  =  [Cc]M  +  mM 

where  {x^}  =  {x*},  {uc}  =  {Vc}  =  {ye,},  and 


[^] 

= 

1  +  [5e2 

ll/,li<J„  1 

[5c] 

=  [5e, 

1  +  [5e2 

[a] 

=  [a, 

J  +  ] 

[Dc] 

a. 

=  ID.JIC,,  ] 

~  ['^e22][-^e2i] 

^€22 

=  [/- 

■5e22  ^2,  “ 

-1 

(4.24) 


(4.25) 


The  well-posedness  of  the  connections  in  Eq.  (4.23)  guarantees  that  [/  -  12]“^  exists. 

The  model  structure  and  its  parameters/matrixes  in  Eqs.  (4.24)  and  (4.25)  axe  indepen¬ 
dent  of  the  structural  and  aerodynamic  variables  and  thus  are  unchanged  in  the  design  or 
optimization  process  of  the  latter 
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4.4  The  ASE  Model 


The  ASE  model  is  obtained  by  connecting  the  AE  model  of  Eqs.  (4.14)  and  (4.15)  with  the 
control  system  of  Eq.  (4.24)  through  fixed  and  variable  gain  connections,  as  shown  schemat¬ 
ically  in  Figure  4.2.  As  in  subsection  4.3.2,  an  input  and  output  of  one  fixed  connection 
cannot  appear  in  another  connection. 


AE  Plant 


I - , 

Control  System 

Figure  4.2:  The  ASE  interconnection  model.  Thick  lines  represent  fixed  connections. 
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To  perform  these  interconnections,  the  system  equations  (4.14),  (4.15)  and  (4.24)  are 
assembled  and  their  inputs  and  outputs  are  re-arranged  as  follows: 


/  1  _  ^  \  4-  "^P*  ^  ®  ^P2 

\  Sc  j  [  0  >lc  J  \  Sc  /  [  0  0  .Bcj  Bcj  J  Uci 

Uc2 

{vpi  ]  r  cu  0  1  r  0  0  0  0  ]  f  Up, 

yp2  ^  _  ^P2  0  /®p1_l  ^  ®  ®  ®  ;  “P2 

Vci  I  0  Cc,  \  Sc  /  0  0  Den  I 

yC2  ,  .  0  Cc2  .  _  0  0  Dc2j  ■^Cjj  .  I  ^C2 

where  inputs  and  outputs  with  the  subscript  1  are  used  in  connection  through  the  vari¬ 
able  gain  matrix,  while  the  inputs  and  outputs  with  the  subscript  2  are  used  in  the  fixed 
connections,  performed  first.  The  system  with  fixed  connections  closed  but  variable  gain 
connections  open  will  be  referred  to  as  the  gain-open  system. 

The  fixed  connections  are  specified  by 


f  ^P2  1  [  0  4c  1  f  J/P2  1 

(“02/  -fep  ®  .  1  2/c2  / 


(4.27) 


where  [/pc]  and  [/cp]  are  of  a  similar  structure  as  [/2]  of  Eq.  (4.23).  Eq.  (4.27)  expresses  fixed 
connections  between  the  AE  and  control  model.  Such  connections  within  the  control  model 


were  already  implemented  in  Eq.  (4.24),  while  fixed  connections  within  the  AE  model  are 
not  logical  and  thus  forbidden.  (Any  connections  from  the  sensors  to  actuators  are  assumed 
to  be  part  of  the  control  system  and  thus  performed  through  the  elements  of  the  latter.) 
Substituting  Eq.  (4.27)  in  Eq.  (4.26)  yields  the  gain-open  vehicle  equations 


{®„}  =  [A«]{x„}  -I-  [D„]{u„} 
{j/v}  ~  "I" 


where 
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and 


[^]  = 
[a]  = 


i4p  “h  Bp2lpcDc22^Cp(^p2 

Bp2  IpcGc2 

{B.]  = 

Bp2  IpcDc2i 

Bc2lcpGp2 

A, 

0 

Be, 

Gp, 

Dci2  ^CpOp2 

1 - 1 

m  = 

o  o 

0 

Ben 

(4.29) 


It  can  be  observed  that  in  Eq.  (4.29),  direct  input-to-output  feed-through,  i.e.,  non-zero 
entries  in  [D„],  can  exist  only  between  Uc^  and  j/ci ,  namely  [Pen]  coimecting  between  internal 
control  elements.  This  result  is  due  to  the  zero  feed-through  of  the  AE  system  discussed  in 
section  4.2,  Eq.  (4.15). 

The  final  ASE  variable  gain  loop  is  closed  by  relating  the  input  vector  {t^v}  to  the  output 
vector  {y„}  via  a  gain  matrix  [(?«], 


{u„}  =  [(?v]{j/v} 


(4.30) 


or  more  specifically 


I  I  =  f  2"^  ]  I  (4.31) 

I  i  L  J  I  i 

Gpp  represents  a  direct  connection  through  a  variable  gain  from  AE  system  sensors  and 
actuators,  while  Gec  are  control  system  internal  variable  gains.  Gep  and  Gpe  are  variable 
gains  specified  on  the  inputs  and  outputs,  respectively,  of  the  control  system,  i.e.,  gains  on 
the  sensors  and  to  the  actuators. 

Substituting  Elq.  (4.30)  into  Eq.  (4.29)  yields  the  closed-loop  ASE  equations  of  motion 


(4.32) 


where 


[A,]  =  [A>]  +  [B„][G,]  [I  -  D,G,r^  [C7„]  (4.33) 

The  inverse  in  Eq.  (4.33)  exists  assuming  the  feedback  given  by  Eq.  (4.30)  is  well  posed. 

The  special  structure  of  the  coefficient  matrices  in  Eq.  (4.29)  is  used  to  reduce  the  com¬ 
putational  load  in  constructing  [Ay].  In  addition,  in  AE  optimization,  only  elements  of  the 
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AE  model,  i.e.,  matrices  with  the  subscript  p  in  Eq.  (4.29),  have  to  be  recomputed,  while 
the  rest  of  the  terms  in  the  [A„],  [B„],  \G^  and  [£)„]  matrices  are  unchanged.  This  can  also 
reduce  the  computational  load  in  the  optimization  process. 


4.5  System  Matrix  Derivatives 


The  derivative  of  \Ay\  in  Eq.  (4.33)  with  respect  to  a  structural  design  variable  Vi  is 


(4.34) 


The  only  terms  in  [>1„]  which  are  functions  of  the  structural  design  variables  are  in 
which  combines  the  second  row  partitions  of  [Aae]  and  [Ba^  of  Eq.  (4.5),  namely 

[A^j  =  -[M]-^[i(2)]  (435) 

where 

=  [  Khh  +  Bhh  +  ^Ahhi  qDh  qAhco  ^A/,ci  ] 

The  derivative  of  with  respect  to  Vi  is 

([DGMV]i[A(2)]  +  [  [DGKV]i  0  0  0  0  [DGMVC]i  ])  (4.36) 

where  [DGMVji,  [DGKV].-  and  [DGMVC]i  are  defined  in  Eqs.  (2.7),  (2.9),  and  (2.13),  and 
where  use  is  made  of  the  differentiation  formula 


(4,37) 


dvi  ^  dvi 

The  derivative  of  [A^]  with  respect  to  Vi  is  all  zero  except  for  the  first  2n/,  +  Ua  columns  in 
rows  rzfc  +  1  to  2nh  which  are  given  by  Eq.  (4.36) . 

The  derivative  of  [Cu]  with  respect  to  Vi  is  zero  in  the  cases  of  displacement  and  velocity 
measurements.  In  the  case  of  acceleration  measurements,  Elq.  (4.8)  yields 


dK] 

dvi 


I  0 

0  i?ci2  lep 


[<^y] 


(4.38) 
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4.6  Aeroelastic  System  Gain  and  Its  Sensitivities 

The  aeroelastic  system  gain  is  useful  in  specifying  the  closed  loop  ASE  system  performance. 
In  particular,  the  zero-frequency  (DC)  gain  of  the  AE  system  excluding  the  rigid  body 
modes,  often  referred  to  as  the  actuator  effectiveness,  can  be  used  to  specify  low-frequency 
performance.  Using  Eqs.  (4.14)  and  (4.15),  the  actuator  effectiveness  is  computed  by 

[<3p(0)l  =  -lC,)K|-‘[Bpl  (4.39) 


Its  sensitivity  with  respect  to  a  structural  design  variable  Vi  is 


9[G^(0)l 

dvi 


\C,]\A,] 


_.a[A,i  d[c,\ 


dvi 


dvi 


(4.40) 


where  the  sensitivity  matrices  9[Ap]/5u,-  and  5[C'p]/9uj  were  discussed  in  the  previous  section. 
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Chapter  5 

Flutter  Analysis 


5.1  Flutter  Constraints 


The  flutter  analysis  is  based  on  the  TVoot  eigenvalues  A{,  whose  Im(A/)  >  0,  of  the  closed-loop 
systm  matrix  [Av]  of  Eq.  (4.33).  Without  a  control  system,  [A]  =  [Aie]  where  [Aoe]  is  the 
open-loop  system  matrix  defined  in  Eq.  (4.5). 

The  flutter  constraints  are  those  specified  in  the  standard  FLUTTER  discipline  in  AS¬ 
TROS,  namely 


^  - IjRBQ  ^  ^  i  =  l,2,...,nu 

^  GFACT  -  Z=  1,2,..., 


(5.1) 


where 

Re(A,) 

^  Im(Ai) 

at  the  yth  user-defined  velocity  associated  with  a  Mach-density  pair.  The  derivative  of  the 
flutter  constraint  with  respect  to  a  design  variable  Vi  is 


dg  ^  1  d'iji  ^  DRRV,  -  7,7DIRV^  ,  . 

dvi  GFACT  dvi  Im(Aj)GFACT 

where  DRRVi  and  DIRVi  are  the  derivatives  of  the  real  and  imaginary  parts  of  d^ijdvi. 
Derivatives  are  calculated  for  a  small  portion  of  the  extracted  eigenvalues  which  axe  the  most 
critical.  The  computation  of  dXifdvi  requires  the  extraction  of  the  eigenvectors  associated 
with  Aj. 
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5.2  Eigenvalue  Sensitivity 

The  complex  column  (left)  eigenvector  {X/}  and  row  (right)  eigenvector  associated 


with  Xi  satisfy 


([i.)  -  A,[/l)  {:?,}  =  {0} 


((A,(„)1-A,[/1)={oF  (6.4) 

To  compute  the  flutter  column  eigenvector  efficiently,  Eq.  (5.3)  is  pautitioned  into 

/[  0  ^  0  1  \  f  ' 

I  -^21  A22  A2i  —  Aj[/]  I  ■!  X2  y  =  {0}  (5.5) 

\L  -^31  A32  Az3  J  /[Xs 

where  subscript  1  relates  to  the  structural  displacement  states  {^},  2  relates  to  the  velocity 
states  {^},  and  3  relates  to  the  remaining  states.  The  first  row  paxtition  in  Eq.  (5.5)  yields 

{X.}  =  (5.6) 

which,  when  substituted  in  the  remaining  parts  of  Elq.  (5.5),  yields 


A22  +  A2Z 

A32  +  ^-^31  -^33 


=  {0} 


The  first  term  in  {Xj}  is  set  to 


X2,  =  (1.0, 0.0) 


The  other  values  of  {^2}  and  {X3},  combined  in  {X},  are  found  by  solving 

[lA]  -  A,[/])  {X}  =  -{6}  (6.9) 

where  [.A]  is  the  partitioned  matrix  in  the  right  side  of  Eq.  (5.7),  without  the  first  row  and 
first  column,  and  {6}  is  the  first  column  of  the  partitioned  matrix  (without  the  first  term). 
Eq..  (5.9)  is  solved  by  first  decomposing  [A]  —  Aj[/]  into 


[A]  -  A,[/]  =  [L][U] 


(5.10) 
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where  [L]  and  [C/^]  are  lower  and  upper  triangular  matrices.  Forward-backward  substitution 
is  then  performed  to  yield  {X},  which  combines  with  Eqs.  (5.6)  and  (5.8)  to  the  column 
eigenvector  {Xj}. 

For  the  computation  of  the  flutter  row  eigenvector,  Eq.  (5.4)  is  partitioned  into 


■  0  I  _0 
-^21  -^22  -423 
.  -Asi  .432  i433 

The  first  column  partition  of  Eq.  (5.11)  yields 


=  {0}^ 


(5.11) 


{X.}  =  i  +  [X3jj’'{X3}) 


(6.12) 


which,  when  substituted  in  the  remaining  parts  of  Eq.  (5.11),  yields 


{*)'( 


■422  +  -^-421  A23 


A32  +  3:7 -431  A33 


-  A,[/]  =  {0} 


(5.13) 


Similarly  to  the  column-eigenvector  solution,  The  first  term  in  {A’2}  is  set  to 


X2,  =(1.0, 0.0) 


(5.14) 


and  the  other  values  of  and  {{-<^3},  combined  in  can  now  be  found  by  solving 

{Xf  ([i]  -  A,(/])  =  -{if  (5.16) 

where  the  left-side  coefficient  matrix  is  identical  to  that  of  Eq.  (5.9)  but  {6}  now  contains 
the  first  row  of  the  partitioned  matrix,  without  the  first  term.  In  order  to  solve  Eq.  (5.15) 
with  the  same  decomposition  used  for  solving  Eq.  (5.9),the  problem  is  posed  as  the  conjugate 
transpose  of  Eq.  (5.15), 

(l-4-)-A?[/j)’'{X-}  =  -{i'}  (5.16) 

from  which  and  then  the  other  parts  of  in  E^qs.  (5.12)  and  (5.14)  are  recovered. 
The  differentiation  of  Eq.  (5.3)  with  respect  to  a  design  variable  Vi  yields 

^  =  W  (5.17) 
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While  the  optimization  flutter  constraints  are  based  on  the  system  eigenvalues  at  user- 
defined  flight  conditions,  the  final  analysis  may  include  the  comutation  of  the  flutter  velocity 
at  which  a  root  branch  crosses  to  the  right  side  of  the  Laplace  domain.  In  addition  to 
the  n„ef  major  velocity  values  (V)  defined  by  the  user  for  each  Mach-density  case,  the  user 
defines  intermediate  V  values  by  specifying  the  number  of  equal  intervals  {njei,)  between  two 
consequtive  major  V  values.  Startig  with  the  third  V  value  in  each  Mach-density  case,  the 
flutter  velocity  is  calculated  when  a  real  part  of  a  complex  eigencvalue  becomes  possitive. 
The  total  number  of  V  values  for  which  roots  are  calculated  in  the  final  analysis  is  nvxndv+1. 
The  kth.  intermediate  V  between  the  major  Vm  and  Kn+i  is 

Vk  =  qm+  —  -  1)  (5.19) 

Udev 

The  flutter  velocity  and  frequency,  Vf  and  Wf  are  calculated  by  quadratic  interpolation  of 
the  first  eigenvalue  that  crosses  to  the  right  side,  and  the  associated  ones  in  the  two  previous 
V  points.  The  eigenvalues  that  correspond  to  the  same  ”  branch”  may  not  appear  in  the 
same  location  in  the  vectors  of  eigenvalues  associated  with  different  V’s.  If,  for  example, 
\ki  is  the  ith  eigenvalue  in  the  vector  of  eigenvalues  {A}jfe  associate  with  V*,  may  not 

relate  to  the  same  aeroelastic  branch.  It  is  assumed,  however,  that  the  V  increments  are 
small  such  that,  the  same-branch  eigenvalue  is  the  closest  one  among  A(ji,_i)^i_2)  to  A(fc_i)jj.j.2). 
The  flutter  boundary  routine  finds  the  first  eigenvalue  in  each  branch  (if  any),  Zfc  +  that 
has  positive  real  part,  and  the  two  previous  same-branch  eigenvectors,  Zk-i  -f-  Wk-i  and 
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Zk-2  +  io^k-2-  The  two  quadratic  interpolation  fornaula  are 


V  =  Vf  +  Agz  +  Bgz^  (5.20) 

uj  =  cjf  +  A^z  4- 


yield  the  Vandermonde  equations 


1 

1 

_  1 

and 

■  1 
1 

_  1 

which  yield 


Zk-2 

Zk-1 

r2  1 

^k-2 

z2 

^k-l 

fl'l 

^  _  i 

!  ) 

Zk 

^k  J 

V4  J 

Zk-2 

^k-2 

{  ^k-2 

Zk-1 

4-1 

A^ 

>  =  <  <^k-l 

Zk 

4  . 

\  J 

1  . 

(5.21) 


(5.22) 


Vk-2Zk-lZk(zk  -  Zk-l)  -  Vk-lZk-2Zk{zk  -  Zk-2)  +  VkZk-2Zk-l{zk-l  -  Zk-2) 
i^k  ^k—l){zk  Zk—2)(^Zk—l  Zk^2) 


(5.23) 


and 


OJf 


^k-2Zk-\Zk[zk  —  Zk^i)  -  U)k-lZk-2Zk{Zk  -  Zk-2)  +  UJkZk-2Zk-l{zk-l  -  Zk-2) 
{Zk  -  Zk-l){zk  -  Zk-2){Zk-lZk-2) 


(5.24) 
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Chapter  6 

Control  Stability  Margins 


This  Chapter  discusses  the  controller  stability  margins  of  the  closed  loop  ASE  system  de¬ 
scribed  in  Chapter  4.  The  analysis  assumes  variations  in  the  gain  matrix  [(?«]  of  Eq.  (4.30) 
only,  applying  both  Single-Input-Single-Output  (SISO)  and  Multi-Input-Multi-Output 
(MIMO)  techniques.  The  SISO  stability  margin  discussions  address  the  single  gain  un¬ 
certainty  cases.  This  includes  descriptions  of  methods  to  calculate  the  appropriate  SISO 
transfer  functions  from  the  MIMO  control  system,  as  well  as  methods  to  calculate  SISO 
gain  and  phase  margins  and  their  sensitivities  with  respect  to  variations  in  structural  design 
variables.  The  SISO  techniques  are  then  further  extended  to  treat  MIMO  problems,  where 
several  controller  gains  can  vary  simultaneously. 

6.1  SISO  Stability  Margins 

The  SISO  analysis  addresses  the  case  of  uncertainties  or  possible  variations  in  one  element 
of  the  controller  gain  matrix  [(j„],  e.g.,  ,  which  takes  the  form 

(6.1) 

Here,  Gvij  is  the  nominal  values  of  ^ ,  while  g-^  and  ^  are  its  gain  and  phase  variations, 
respectively,  which  determine  the  gain  and  phase  margins.  Their  nominal  values  are  =  1 
and  ^  =  0.  Assuming  that  nominally  the  closed  loop  ASE  system  is  asymptotically  stable. 
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i.e.,  [^„]  of  Eq.  (4.33)  is  Hurwitz  with  the  nominal  ^  ,  SISO  stability  margins  determine 
the  smallest  variations  in  g^  .  and  ^  that  cause  instability  of  the  closed  loop  ASE  system. 
To  compute  these  margins,  appropriate  SISO  transfer  functions  have  to  be  constructed  from 
the  (generally)  MIMO  control  system.  Consequently,  the  SISO  gain  and  phase  margins  and 
their  sensitivities  to  structural  design  variables  can  be  computed. 

6.1.1  SISO  Transfer  Functions  from  a  MIMO  Control  System 

Computation  of  the  SISO  stability  margins  against  variations  in  discussed  above  can 
be  efficiently  performed  by  opening  the  i-j  loop  of  the  MIMO  control  loop  in  Eq.  (4.30)  at 
either  input  or  output  side  of  the  gain  element  G„.  ^ ,  resulting  in  a  new  single  variable  input 
and  a  new  single  variable  output  to  the  system.  The  SISO  stability  margins  can  then  be 
computed  from  the  open  loop  SISO  transfer  function  between  this  new  input  and  the  new 
output. 


Figure  6.1:  SISO  stability  margins:  control  loop  layout  for  the  required 
SISO  transfer  function. 


This  new  SISO  transfer  fimction  is  computed  using  the  “open-closed”  loop  system  de¬ 
picted  schematically  in  Figure  6.1,  where  the  elements  of  the  new  gain  matrix  [G(’’^)=°]  are 
equal  to  the  elements  of  the  nominal  gain  matrix  [Gv],  except  for  the  {i,j)  element  which  is 
set  to  zero  and  {e,}  and  {ej}  are  the  i-th  and  the  j-th.  unit  vectors,  rj  is  the  new  “external” 
scalar  input  and  yj  is  the  new  input  of  this  open-closed  system. 
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The  MIMO  loop  in  Fig.  6.1  is  defined  by 


{“•}  =  [gI’’’*  “]  {s/«}  +  {ei}n 


(6.2) 


Substituting  the  controller  Eq.  (6.2)  into  the  system  Eq.  (4.28),  and  using  the  relation 


a;  =  }’■{».} 


(6.3) 


yields  the  new  open-closed  loop  system  state-space  equations 

■{•Coc}  —  [.^ocj'f^oc}  d” 


where 


Vj  —  “I" 

14J  =  [A.]  +  [b.]  [C,] 

{W  =  IB.1{^}+[b.]{b.,} 
hi]  =  G„,,{e,F[C.] 

[b.]  =  |B.]  [(?!)•>)=»] 

6,]‘'  lD.|{ei} 


{d..}  = 
[B.]  = 


The  SISO  tramsfer  function  between  and  j/y  is 


Vjis)  =  Toc^A^)  ri(s) 

ToCj^('S)  =  [■'^oc]]  {^OCi}  d"  d< 


(6.4) 


(6.5) 


(6.6) 

(6.7) 


This  transfer  function  is  used  next  to  compute  the  gain  eind  phase  margins  of  the  control 
system  gain,  and  their  corresponding  sensitivities  to  variations  in  the  structural  parameters 
of  the  ASE  system. 
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6.1.2  SISO  Gain  Margins 

Upper  and  lower  gain  margins  are  defined  for  a  stable  closed  loop  ASE  system.  The  upper 
gain  margin  indicates  the  lowest  relative  increase  in  the  loop  gain  ^  for  which  the  closed 
loop  system  is  imstable.  Using  the  notation  of  Eq.  (6.1),  the  upper  gain  margin  is  the  smallest 
9x^  >  1  that  causes  closed  loop  instability.  Similarly,  the  lower  gain  margin  indicates  the 
lowest  relative  decrease  in  i.e.,  the  largest  0  <  <  1,  that  leads  to  an  unstable  closed 

loop  system.  According  to  the  Nyquist  stability  criterion,  gain  margins  are  the  reciprocals 
of  the  open  loop  transfer  function  gain,  evaluated  at  the  phase  crossover  frequencies  defined 
as  the  frequencies  for  which  the  transfer  function  phase  equals  to  0°.  The  0°  (and  not 
the  classical  -180°)  phase  crossing  is  used  because  of  the  positive  feedback  assumption  in 
Eq.  (4.30).  Among  all  the  phase  crossover  points,  the  two  closest  to  +1  in  the  complex 
Nyquist  plane  determine  the  upper  and  lower  gain  margins. 

In  the  ASE  context,  the  SISO  transfer  function  of  interest  is  .(«)  defined  in  Eq.  (6.7). 
To  compute  the  gain  margins,  it  is  necessary  to  determine  the  phase  crossover  frequencies  uipco 
for  which  (ja;pco)  is  real  and  positive.  For  that,  the  matrices  [2/]  and  [N]  are  constructed 
as  follows 

[Aoc]  {^oci}  r  ^  - 

-[coc,]  1^]=  ;  (0.8) 

.  N  -{b^f  0  J  [0 . 0  . 

where  x  n„  is  the  dimension  of  [Ax:]-  The  phase  crossover  frequencies  are  the  imaginary 
parts  of  the  purely  imaginary  generalized  eigenvalues  of  [L]  and  [AT],  i.e. 

[L]  {xfc}  =  Afc  [N]  {sjk} ,  A:  =  1, . . . ,  (6.9) 

and  {®fc}  are  the  corresponding  generalized  right  eigenvectors.  Among  all  the  purely  imag¬ 
inary  Afc-s  only  the  ones  with  real  and  positive  Tocj  -{)^k)  relevant  to  finding  upper  and 
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lower  gain  margins,  leading  to  the  set 


a> 


<  Imag(Ak) 


Real  =  0,  Imag  ^Ak)  >  0, 

Real  (Toq^CAk))  >  0,  Imag  (Tocj^(Ak))  =  0 


By  definition,  the  upper  and  lower  gain  margins  are  determined  as 


(6.10) 


i -  (6.11) 

Max  |o  <  <  l} 

GMi  =  - - -  (6.12) 

Min  {l  <  Tocj,i{ju>)} 


The  expressions  above  yield  the  phase  crossover  frequencies,  and  which  correspond 
to  the  upper  gain  margin,  (?Mu,  and  lower  gain  margin,  GMi,  respectively.  Therefore, 


and  in  decibels. 


GMu,i  = 


2oCj,< 


u>I 


C7Il4,,[dB]  =  201ogGMu,i=  -201ogToq,(H| 


(6.13) 


(6.14) 


Note:  Some  systems  may  have  no  upper  and/or  lower  gain  margins  (often  set  to  ±oo  [dB]). 
In  these  cases,  the  gain  margins  and  their  corresponding  sensitivities,  discussed  below,  are 
not  calculated. 


6.1.3  SISO  Phase  Margins 

Upper  and  lower  phase  margins  are  defined  for  a  stable  closed  loop  ASE  system.  The  upper 
phcise  margin  indicates  the  lowest  additional  loop  phase  lag  for  which  the  closed  loop  system 
is  unstable.  Using  the  notation  of  Eq.  (6.1),  the  upper  phase  margin  is  minus  the  largest 
<f>ij  <  0  that  causes  closed  loop  instability.  Similarly,  the  lower  phase  margin  indicates  the 
lowest  additional  loop  phase  lead,  i.e.,  minus  the  smallest  <t>ij  >  0,  that  leads  to  an  unstable 
closed  loop  system.  Based  on  the  Nyquist  stability  criterion  for  positive  feedback  systems. 


55 


phase  margins  are  evaluated  at  the  gain  crossover  frequencies  (defined  as  the  frequencies 
for  which  the  transfer  function  gain  is  unity)  and  equal  to  the  transfer  function  phase  at 
these  frequencies.  Among  all  the  gain  crossover  points,  the  two  closest  to  +1  in  the  complex 
Nyquist  plane  are  used  to  determine  the  upper  and  lower  phase  margins.  By  convention,  the 
upper  phase  margin  is  positive  and  the  lower  phase  margin  is  negative. 

In  the  ASE  context,  of  Eq.  (6.7)  is  used  as  the  open  loop  transfer  function  for 

phase  margin  calculations.  The  gain  crossover  frequencies  ojgco  Q-re  found  by  computing  the 
eigenvalues  of  the  Hamiltonian  matrix  [/f]  defined  by 


[H]^ 


[Aqc]  “h  ^OCj 


1  dn. 


{V,}{6„,}  /4 


T 

~  [cocj  ]  [^ocj]  /dj,i 


1  ’’ 


[■Aoc]  +  {6oci }  [Cj\ 

Hi 


(6.15) 


where  dj^i  —  1  —  Among  these  eigenvalues  A*,  i  =  k,. . . ,  ,  only  purely  imaginary 

ones  for  which  rocj  j(A*.)j  =  1  are  used  to  determine  the  upper  and  lower  phase  margins, 
leading  to  the  set 


a;  :=  I  Imag(Ak)  |  Real  (At)  =  0,  Imag  (A^)  >  0,  |Tocj,i(Ak)|  =  l|  (6.16) 
The  upper  and  lower  phase  margins  axe  determined  as 


=  Min  |o°  <  /  Tocj.i  (jfa^)  <  180°}  (6.17) 

PMi  =  Max  {-180°  <  <  0°)  (6.18) 

where  is  manipulated  to  be  in  the  range  -180°  <  /7’oc,,<(j^)  <  180°.  Note  that 

the  above  expressions  yield  the  conventionally  assumed  >  0  and  PMi  <  0.  These 

expressions  also  yield  the  gain  crossover  frequencies,  and  which  correspond  to  the 
upper  phase  margin,  PM„,  and  lower  phase  margin,  PM/,  respectively.  Therefore, 


PM«,/  =  /Toc,-,(ja;) 


(6.19) 
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PMu,i  =  tan' 


-1  /Imag[Toq,,Oa;)] 


(6.20) 


VReal[To,,(ja,)]; 

Wgco 

Note:  Some  systems  may  have  no  upper  and/or  lower  phase  margins  (often  set  to  ±oo[°]). 
In  these  cases,  the  phase  margins  and  their  corresponding  sensitivities,  discussed  below,  are 
not  calculated. 

6.1.4  Sensitivities  of  SISO  Gain  Margins 

In  this  section  the  sensitivities  of  the  upper  and  lower  gain  margins  with  respect  to  a  design 
variable  are  derived.  The  gain  margins  axe  computed  for  the  controller  gain  variations, 
while  the  design  variables  include  structural  variables  and  the  controller  variable  gain  matrix 
[G„]  entries,  all  grouped  in  the  vector  {v}. 

The  sensitivity  of  the  gain  margins  expressed  in  decibels  are  given  by 


dGM^jjdB]  _  20/ In  10  dGMu,i 
dvi  GMu,i  dvi 

where  dGM^^ijdvi  is  the  regular  (not  in  decibel)  gain  margin  sensitivity. 
The  gain  margin  sensitivities  with  respect  to  ,  are  given  simply  by 

dGM^,i  __  GM^,i 

dG.,,  a.. 


(6.21) 


(6.22) 


or  in  decibels 

aGM„, ,[dB]_  20/ In  10  .  . 

a,,  ^ 

To  compute  the  gain  margin  sensitivities  with  respect  to  the  other  design  variables 
the  matrix  [4GAfu.i]  is  constructed 


(6.24) 
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Note  that  is  an  eigenvalue  of  j^AGrM„,i]-  The  gain  margin  sensitivities  are  given  by 

w'l 


dGM^j 

dvi 


dvj 

{n"  {jf } 


R«1  - - 

{yy  {xy 


T 

where  {X}  and  {F}  are  right  and  left  eigenvectors  of  given  by 

{S'}”  A™.,.]  =  K^'Myf 


(6.25) 


(6.26) 


(6.27) 


f^OCj 


(6-2«) 

The  sensitivity  dAGM^^Jdvi  is  constructed  using  the  results  of  Section  4.5. 

6.1.5  Sensitivities  of  SISO  Phase  Margins 

Upper  and  lower  phase  margin  sensitivities  are  computed  with  respect  to  the  design  variable 
vector  {u}  defined  in  the  previous  section.  To  compute  these  sensitivities,  the  matrix  \Apm^  ,] 
is  constructed 


=  [Aoc]  + 


gjPAfu,i[rad]  _  ^ 


{^oci}  [cocj] 


(6.29) 


where  PM„,/[rad]  is  the  phase  margin  expressed  in  radians.  Note  that  is  complex 

and  is  its  eigenvalue. 
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The  phase  margin  sensitivities  are  given  by 


Real 


dPM^,i  _ 

dvi 


oTw 


Real 


aPM„,,[rad] 


{^} 


_  m 
\ 


{Y}  {X} 


where  {X}  and  {1^}^  axe  right  and  left  eigenvectors  of  [j4pAf„  ,]  given  by 


{F}  =  J(^gco{Y} 


and 


dApM^, 


je 


7PjWu^i[rad] 


{^oci}  [^ocj 


aPMu, /[rad]  y 

Again,  the  sensitivity  dApM^Jdvi  is  constructed  using  the  results  of  Section  4.5. 


(6.30) 


(6.31) 

(6.32) 

(6.33) 


6.2  MIMO  Stability  Margins  and  Sensitivities 


Following  Ref.  14,  the  MIMO  stability  margins  are  defined  by  introducing  ASE  system 
uncertainties  at  the  plant  input  or  output  as  described  in  Figure  6.2. 

The  input  and  output  MIMO  stability  margins  are  defined  with  respect  to  uncertainties 
(variations)  in  the  matrices  [2/^]  and  [Lo],  which  are  assumed  to  be  diagonal,  i.e. 

[LJ  =  diag  [s*e^*i] ,  l  =  l,2,...,Ni  (6.34) 

[tj  =  diag[sje^2],  fc=l,2 . N,  (6.35) 

where  Ni  and  No  are  the  number  of  inputs  and  outputs,  respectively,  to  the  ASE  model. 
At  the  nominal  condition,  g]  =  =  1  and  <l>]  =  <j>%  =  0  ^  I,  k,  i.e.,  [L/j  = 
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Figure  6.2:  MIMO  stability  margin  evaluation  setup. 


[Xto]  —  Sufficient  conditions  for  a  guaranteed  gain  margin  g*  and  phase  margin  at 

any  input  to  the  gain-open  ASE  plant  are 


[(1  -  1/^)2  -h  2(1  -  cos0')/p‘]'^'  <  5l{/  -  (6.36) 

[(1  -  +  2p’(l  -  cos.^')]'^'  <o:{l-  (6.37) 


where 


[«.(»)]  =  [ajK-A.i-‘[Bj+iD.] 


and  £{•}  is  the  minimum  singular  value  of  the  argument  matrix  evaluated  at  s  =  ja;  V  a;  >  0. 
Sufficient  conditions  for  output  gain  and  phase  margins  are  similar  to  the  inequalities  of 
Eqs.  (6.36)  and  (6.37)  except  of  replacing  the  [G^Pv]  terms  with  [P„G„]. 

The  inequalities  of  Eqs.  (6.36)  and  (6.37)  on  the  ASE  system  minimum  singular  values 
can  serve  as  MIMO  stability  margins  in  the  multidisciplinary  optimization  process.  Their 
sensitivities  with  respect  to  a  design  parameter  p  are  computed  analytically  from  the  singular 
value  derivative  expression 


da{H} 

dp 


=  Real 


(6.38) 
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where  u  and  v  are,  respectively,  the  right  and  left  normalized  singular  vectors  of  the  matrix 
[H],  corresponding  to  the  minimum  singular  value  £.{H}.  In  the  ASE  MIMO  case, 

[fl]  =  [/]  -  [G.]  [tC.Ks/  -  +  P.l]  (6.39) 

and  the  evaluation  of  the  derivative  da{H}jdp  facilitates  the  sensitivities  of  [Ai,],  [B„],  [(7v] 
and  [D„]  which  were  discussed  earlier  in  the  SISO  margin  section. 
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Chapter  7 

Continuous  Gust  Response 


The  time-domain  ASE  model  for  continuous  gust  response  analysis  is  is  based  on  the  ASE 
model  of  Chapter  4,  augmented  by  gust  states  and  excited  by  a  gust  velocity  input.  The 
formulation  is  based  on  the  gust  modeling  of  Ref.  5.  It  assumes  a  single  gust  pattern,  either 
vertical  or  lateral,  represented  by  the  spectral  properties  of  the  gust  velocity  amplitude. 
The  gusts  are  assumed  to  be  spanwise  uniform.  Kinematic  load  modes  and  the  associated 
aerodynamic  and  inertial  coefficient  matrices  facilitate  the  use  of  section  loads  as  design 
constraints. 

7.1  Equation  of  Motion 

The  second-order  time-domain  equation  of  motion  of  the  structure  was  defined  in  Eq.  (2.23), 
with  the  generalized  external  forces  represented  by  {/),(<)}.  The  external  forces  for  stability 
analysis  were  expressed  in  Section  4,  where  they  included  aerodynamic  forces  due  to  struc¬ 
tural  dynamics,  and  inertial  and  aerodynamic  forces  due  to  control-surface  motion.  These 
loads  can  be  considered  as  part  of  the  internal  dynamics  of  the  closed-loop  aeroservoelastic 
system.  Gust  response  analysis  requires  the  addition  of  gust  input  loads.  When  output 
section  loads  are  required,  generalized  aerodynamic  loads  associated  with  kinematic  load 
modes  are  required  as  well. 
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The  generalized  aerodynamic  forces  are  defined  in  the  frequency  domain  in  Eq.  (2.38). 
The  aerodynamic  force  coeflhcient  (AFC)  matrices,  including  the  gust  column  ‘•■re 

first  calculated  at  several  user-defined  tabulated  reduced  frequency  values,  kt.  The  tabulated 
matrices  axe  used  for  approximating  the  AFC  matrix  as  a  rational  function  of  k  in  the  entire 
frequency  domain,  as  described  in  Chapter  3.  An  expansion  to  the  entire  Laplace  domain 
is  performed  by  replacing  ik  in  the  rational  expression  by  the  non-dimensional  Laplace 
variable  p  =  shjVy  which  yields  Eq.  (3.2).  The  substitution  of  p  =  s6/V  in  the  expression 
for  [^h(p)]  of  Eq.  (3.2),  and  the  subsequent  construction  of  the  state-space  equations  for 
stability  analysis  were  discussed  in  Chapter  4.  The  expanded  version  of  Eq.  (4.1),  which 
includes  the  gust  related  coefficients,  reads 
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Qh{s) 

Ql{s)  J 


Aho 

Alo 


b 


Ahi 

Ali 


s  + 


y2 


Ah2 

Al2 


Dh 

Dl 


(W»- (71) 


where  the  coefficient  matrices  are  column  partitioned  as 

(n  =  0,l,2). 


Ahn 

Ahhr,  Ahc^  AhOn 

Ain 

_  AlHt,*  Atcn  ALGn 

[E]  =  [Eh  E,  Eg] 


When  the  design  values  are  different  than  the  ones  for  which  the  approximation  matrices 
were  calculated,  {AhOn}  is  updated  by  Pre-multiplying  its  baseline  value  by  and  [ALhn] 
is  updated  by  post-multiplying  it  by  [V>],  similarly  to  the  other  matrix  updated  in  Eq.  (4.2). 

To  avoid  coefficients  associated  with  the  second  time  derivative  of  the  gust  velocity. 
Approximation  constraints  should  be  applied  to  the  gust  columns  to  yield  {A/jGj}  =  0  and 
{A1G3}  —  0-  As  discussed  in  Chapter  3,  the  approximation  matrices  [A;,„],  [D^]  and  [E] 
are  solved  for  while  ignoring  the  section-loads  data.  The  [A^n]  and  [Dl]  matrices  are  then 
solved  for  with  a  fixed  [£]. 

The  augmenting  aerodynamic  state  vector  is  defined  here  by  its  Laplace  transform  as 

=  (w-*  -  jlJJl)'‘  ([i!l.l{«»)}  +  [«.!{«.(■>)}  +  ^  (7.2) 
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Equations  (3.1),  (7.1)  and  (7.2)  add  a  gust  related  term  to  the  generalized  aerodynamic  force 
vector  {Phis)}  of  Eq.  (4.4), 

+  (7.3) 

Equations  (2.23),  (2.28),  (4.4)  and  (7.3)  yield  the  state-space  aeroelastic  equation  of 
motion  for  open-loop  response  analysis 

{®ae}  =  [■Aae]{a5ae}  +  [^ae]{Woe}  +  [Baw\{wG}  (7.4) 

where  {xae},  {uoe},  [-^ae]  a^d  [Bae]  axe  defined  in  Eq.  (4.5)  and 

0  0 
0 

where  [M]  defined  in  Eq.  (4.5). 

Sensor  readings  of  the  aeroelastic  plant  were  defined  in  Eqs.  (4.6)  to  (4.9).  The  only 
effects  of  the  gust  inputs  on  sensor  readings  are  in  the  case  of  accel^'ation  output,  where  the 
expression  for  ya  in  Eq.  (4.8)  is  supplemented  by 

Vac  =  [CaG]{wG}  (7.5) 

where 

[CaG]  =  -[4>y][M]~^  [  ^{A/,Go}  ^{AhGi}  ] 
which  expands  Eq.  (4.9)  to  become 

{j/ae}  =  [C'ae]{®oe}  +  [f^ae]{Wae}  +  {GGa\{wG}  (7.6) 

where  the  only  non-zero  rows  in  [Cgo]  are  those  associated  with  acceleration  signals,  which 
are  taken  from  [C'og]  of  Eq.  (7.5). 
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The  plant  equations  of  motion,  Eqs.  (4.14)  and  (4.15),  become 

{ip}  =  [>lp]{a:p}  +  [Bp]{up}  +  [.Bpu,]{d)Gr} 

{Vp}  =  [<^p]{®p}  +  [C'caji^c} 


where 


The  vehicle  equations  of  motion,  Eq.  (4.28)  become 


{ffiv}  =  [^t;]{a:v}  +  +  [B^w]{wg} 

{Vv}  =  [Gv]{xv}  +  [D„]{u„}  +  [Cq„]{wg} 


where 


[  Bc,Ccg,  ’ 


[Ggv] 


CGa 

Bc\2  ^Cp^G^ 


The  equation  of  motion  of  the  closed-loop  a«roservoelastic  system  is 


{®«}  —  [■'4t,]{®t,}  "I" 


where  [>1„]  is  the  closed-loop  system  matrix  defined  in  Eq.  (4.33)  and 


[B^\  =  [B^]  +  [S„][G„]  [/  -  [CG^ 


Output  parameters  {t/n}  which  are  of  interest  for  dynamic  response  analysis  are  expressed 
in  terms  of  the  plant  states,  similarly  to  the  sensor  readings  of  Eq.  (7.7), 


{Vr}  =  [<^ph]{®p}  +  [GGaR]{u>G} 


(7.10) 


where,  similarly  to  [Cgo]}  [C'g^h]  =  0  except  for  acceleration  response  rows. 


7.2  Gust  Model 


The  description  of  the  gust  mode  and  the  way  gust  column  in  the  generalized  aerodynamic 
matrix  is  calculated,  axe  given  in  the  theoretical  manual  of  the  ZAERO  module.  We  deal 
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in  this  section  with  the  statistical  properties  of  the  gust  velocity  and  their  realization  in 
state-space  modeling. 


It  is  assumed  that  the  gusts  axe  generated  in  a  stationary,  Gaussian,  stationary  random 
process.  A  continuous  atmospheric  is  defined  by  its  root-mean-square  (RMS)  value  and  power 
spectral  density  (PSD)  function.  RMS  gust  values,  <t„q,  for  aircraft  design  are  defined  by 
aviation  regulations.  A  classiceil  formula  for  the  gust  PSD  is  that  of  Dryden, 


2  _  1  +  3(fa;rg)^ 


l  +  (a;T5)^]' 


(7.11) 


where  Tg  =  LgjV  where  Lg  is  the  scale  of  turbulence  typically  2500  ft  for  aircraft  design. 
This  formula  is  consistent  with  the  relation  between  RMS  value  and  PSD  function  in  the 
theory  of  stochastic  processes, 

^hl-oo  ""I  Jo  (7-12) 

The  expression  of  Dryden  appears  some  time  with  a  tt  in  the  denominator,  which  agrees 
with  the  formula  traditionally  used  by  aeroelasticicins, 

(7.13) 


A  more  modern  PSD  formula,  but  somewhat  more  complicated,  is  that  of  Von  Karman, 

\2 


via 


1  +  I  (1.339a>T,)' 

(u;)  —  cr^„Tg 


[l  +  (1.339i.>T,f 


(7.14) 


To  facilitate  the  use  of  modem  algebraic  tools  for  calculating  the  structural  response  to  ' 
continuous  gust  excitation,  we  want  to  describe  the  gust-response  process  by  a  state-space 
model  excited  by  white  noise.  For  this  purpose  we  need  to  define  a  gust  filter  that,  when 
excited  by  white  noise,  produces  the  gust  PSD  function.  The  transfer  function  of  a  filter 
that  produces  Dryden’s  formula  is 


wa{s)  _  '^^3  +  Tj 


(7.16) 


66 


where  w  is  a.  white-noise  parameter  with  =  1.  A  state-space  realization  of  this  filter  would 
result  with  direct  white-noise  excitation  of  the  aircraft  dynamics,  and  hence  non-converged 
RMS  acceleration  response.  To  avoid  that,  we  add  a  low-pass  filter  which  modifies  Tg{s)  to 
be 

W  =  P-16) 

5  i-  a 

A  state-space  realization  of  Tg[s)  is 


{ij  =  [Ag]{Xg}  +  {Bg}w 
{Wg)  =  [Cg]{Xg} 


(7.17) 


where 


■  0 

1 

0  ■ 

f  “  1 

14]  = 

<M 

1 

W  =  -^  0 

0 

0 

—a 

icj  = 


-^/3T-'/=  (l-2^/3)T-W 


Since  Von  Karman’s  PSD,  Eq.  (7.14),  is  a  non-rational  function  of  w,  it  can  not  be 
modeled  exactly  by  a  transfer  function.  Reference  16  suggested  the  3rd  order  filter 

(1  +  2.618Tgs)(l  4-  0.1298x^5) 


1'g  ^  2.083Tgs)(l  -1-  0.823t3s)(1  +  0.0898x^5) 

which  yields  good  fit  in  the  range  0  <  TjW  <  20.  Hohlit^^  suggested 


(7.18) 


Tg  -  ^WQy/^  ^ 

(1  +  2.187r,s)(l  +  0.1833r,s)(l  +  0.021t,s) 

(l  +  1.339T(,s)(l  +  1.118T,s)(l-l-0.1277T,s)(l-h0.0146Tjs)  ’ 

which  is  good  for  0  <  <  200.  Both  approximations  of  Von  Karman’s  PSD  function 

require  additional  low-pass  filter  to  yield  a  state-space  realization  without  an  output  noise 

term,  as  in  Eq.  (7.17).  With  the  low-pass  filter,  Eq.  (7.18)  yields  a  4-state  filter,  and 

Eq.  (7.19)  yields  a  5-state  one. 
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7,3  Gust  Response  Analysis 

Augmentation  of  the  aeroservoelastic  equation  (7.9)  by  the  gust  filter,  Eq.  (7.17),  yields 


{xi}  =  [Ai]{x}  +  {Bi^}w  (7.20) 

where 

{*.}  =  {::}  [^1= 

Gust  response  parameters  expressed  by  Eq.  (7.10)  become 

{yi}  =  (7.21) 

where 

lCiJ  =  [[Cp«]  [Co.h][C,]] 

The  response  paxameters  include  discrete  structural  response  and  integrated  section  loads. 
The  discrete  parameters  are  displacements,  velocities  and  accelerations  that  were  formulated 
in  Eqs.  (4.6)  to  (4.8)  and  (7.5).  Section  loads  can  also  be  expresses  in  the  form  of  Eq.  (7.21). 
Equations  (2.40,  7.1,  7.2,  7.17)  yield  the  time-domain  aerodynamic  loads 

in)  = 

9  [  -^Lho  y-^Lhl  Djj  Alcq  y2  A£,C2  ^  ^A£ff^C'g  j  {s^l} 

-  <l-^[ALh2]{i}  (7.22) 

where  {^}  can  be  expressed  in  terms  of  {xi}  by  using  the  associated  row  partition  in 
Eq.(7.20). 

Rigid-body  displacement  modes  with  no  aerodynamic  stiffness  (in  X,  Y,  Z  and  cause 
zero  roots.  The  rows  and  columns  associated  with  the  deflection  states  of  these  modes  should 
be  removed  to  avoid  singularity  in  the  solution  of  the  Lyapunov  equations  given  below.  We 


[^] 

0 


{Blu,} 


Ba 
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assume  now  that  the  first  state  in  {ici}  after  this  removai  are  0y  in  the  symmetric  case  and 
$z  in  the  anti-symmetric  case.  Other  rigid-body  motion  modes  that  have  a  zero  eigenvalue 
are  the  symmetric  one  with  By  =  Z/V  and  the  antisymmetric  Bz  =  —Y  jV .  Such  singularity 
is  ehminated  by  calculating  the  associated  eigenvector,  creating  a  transformation  matrix  [T] 
in  which  the  first  column  of  a  unit  matrix  is  replaced  by  the  eigenvector,  and  performing  the 
transformation 

{X2}  =  [T]-^{x^}  (7.23) 

which  yields 

{ffi}  = 

{t,}  =  [C7]{x}  (7.24) 


where  {x}  is  {X2}  with  the  first  term  truncated,  [.A]  is  [^2]  =  with  the  first  row 

and  the  first  column  truncated,  [£„,]  is  [B2w]  =  [r][BiJ  with  the  first  row  truncated,  and 
\C]  is  [G^  =  [C'i][r]  with  the  first  column  truncated. 

A  state  covariance  matrix  [X]  is  defined  by  the  expected  value 

[X]  =  E  [{x}{x}’-]  (7.2B) 

When  w  of  Eq.  (7.20)  represents  a  unit-intensity  white-noise  process,  [X]  satisfies  the  Lya¬ 
punov  equation 

M  [A-]  +  [X]lAf  =  (7.26) 

Efficient  solution  is  obtained  via  Schur  decomposition  [A]  =  [i7][!r][t/]^  where  [i7]  is  a  unitary 
matrix,  [U][U]'^  —  [/],  and  [T]  is  an  upper  triangular  matrix.  A  unique  solution  for  [X]  is 
obtained  when  [A]  has  no  roots  with  real  part  equals  zero.  The  output  covariance  matrix  is 
related  to  [X]  by 

m  B  [[C7){x}{x}’-lC'l’']  =  1C1[X][C7F  (7.27) 
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The  diagonal  of  [V]  contains  the  mean-square  output  response 


=  [Cy][X][Cyf  (7.28) 

where  [C^  is  a  row  in  [C\  of  Eq.  (7.24) 

7.4  Sensitivity  Analysis 

The  derivatives  of  gust  response  parameters  with  respect  to  the  design  variables  are  based 
on  the  sensitivity  of  the  [.4i]  and  [(7i]  in  Eq.  (7.20)  and  (7.21).  These  are  transformed  to 
the  derivatives  of  [^4]  and  [C]  in  the  same  way  [Ai]  and  [Ci]  are  transformed  into  [j4]  and 
\G\  in  Section  7.5.  The  derivatives  of  [Aj]  are  based  on  the  derivatives  of  [.A]  discussed  in 
Section  4.5,  the  derivatives  of  which  are  non  zero  only  with  acceleration  sensors.  The 
derivatives  of  [Ag\  and  [Cg]  and  {Bi„}  are  zero.  The  derivatives  of  [Ci]  are  non  zero  only 
when  accelerations  are  involved,  which  happens  in  acceleration  and  section-loads  responses. 
The  differentiation  of  Eq.  (7.28)  with  respect  to  a  design  variable  Vi  yields 

+  2^[c,]|xi!C,r  (r.29) 

The  derivatives  of  [X]  are  obtained  from  the  differentiation  of  Eq.  (7.26),  which  yields  the 
Lyapunov  equation 

which  is  solved  for  d\X\ldvi,  which  is  then  substituted  in  Eq.  (7.26)  for  the  derivatives  of 
the  mean-square  responses.  The  Schur  decomposition  of  [A]  can  be  used  for  all  responses 
and  their  sensitivity  derivatives. 
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